## 3 Enriched Petrov-Galerkin Natural Element Method  
### Enriched NEM Approximation  
The boundary value problem expressed by Equations (1)-(3) is converted to the weak form according to the virtual work principle: Find \(\mathbf{u}(\mathbf{x})\) such that  
\[\int_{\Omega}\varepsilon(\mathbf{v}):\sigma(\mathbf{u})\,d\Omega=\int_{\Gamma_{N}}\mathbf{ \hat{t}}\cdot\mathbf{v}\,ds \tag{22}\]  
for all the test displacements \(\mathbf{v}(\mathbf{x})\) in the Cartesian coordinate system \(\{x,y\}\). Next, a non-convex NEM grid \(\mathfrak{T}_{\textit{NEM}}\), shown in Figure 3a, is constructed using \(N\) nodes and a number of Delaunay triangles. Then, for a given NEM grid, trial and test displacement fields \(\mathbf{u}(\mathbf{x})\) and \(\mathbf{v}(\mathbf{x})\) for PG-NEM approximation are interpolated as  
\[\mathbf{u}_{h}(\mathbf{x})=\sum_{J=1}^{N}\tilde{\mathbf{u}}_{J}\phi_{J}(\mathbf{x})+k_{1}\! \left[ \begin{array}{c}Q_{11}(\mathbf{x})\\ Q_{12}(\mathbf{x})\end{array} \right]+k_{2}\!\left[ \begin{array}{c}Q_{21}(\mathbf{x} )\\ Q_{22}(\mathbf{x})\end{array} \right] \tag{23}\]  
\[\mathbf{v}_{h}(\mathbf{x})=\sum_{I=1}^{N}\tilde{\mathbf{v}}_{I}\Psi_{I}(\mathbf{x})+\lambda_{ 1}\!\left[ \begin{array}{c}Q_{11}(\mathbf{x})\\ Q_{12}(\mathbf{x})\end{array} \right]+\lambda_{2}\!\left[ \begin{array}{c}Q_{21}( \mathbf{x})\\ Q_{22}(\mathbf{x})\end{array} \right] \tag{24}\]  
with Laplace interpolation functions \(\phi_{J}(\mathbf{x})\) represented in Figure 3b. Here, \(\psi_{I}(\mathbf{x})\) are constant-strain FE basis functions, which are defined over Delaunay triangles. In addition, \(k_{1}\) and \(k_{2}\) are two unknown constants associated with a crack, and \(Q_{1\alpha}(\mathbf{x})\) and \(Q_{2\alpha}(\mathbf{x})\) represent the near-tip singular displacement fields given by  
\[Q_{11}(\mathbf{x})=\frac{1}{2\mu}\sqrt{\frac{r}{2\pi}}cos\!\left(\frac{\partial}{ 2}\right)\!\left[\kappa-1+2 sin^{2}\!\left(\frac{\partial}{2}\right)\right] \tag{25}\]  
\[Q_{12}(\mathbf{x})=\frac{1}{2\mu}\sqrt{\frac{r}{2\pi}}sin\!\left(\frac{\partial}{ 2}\right)\!\left[\kappa+1-2 cos^{2}\!\left(\frac{\partial}{2}\right)\right] \tag{26}\]  
\[Q_{21}(\mathbf{x})=\frac{1}{2\mu}\sqrt{\frac{r}{2\pi}}sin\!\left(\frac{\partial}{ 2}\right)\!\left[\kappa+1+2 cos^{2}\!\left(\frac{\partial}{2}\right)\right] \tag{27}\]  
\[Q_{22}(\mathbf{x})=\frac{1}{2\mu}\sqrt{\frac{r}{2\pi}}cos\!\left(\frac{\partial}{ 2}\right)\!\left[\kappa-1-2 sin^{2}\!\left(\frac{\partial}{2}\right)\right] \tag{28}\]Referring to Figure 1, \(r\) is the radial distance from the crack tip, while \(\theta\) is the angle from the \(x^{\prime}-\)axis. Meanwhile, \(\mu\) indicates the shear modulus, and \(\kappa\) is the Kolosov constant given by  
\[\kappa= \begin{cases}3-4\
u&\text{for plane strain}\\ (3-\
u)/(1+\
u)&\text{for plane stress}\end{cases}  \tag{29}\]  
with \(\
u\) being the [PERSON]'s ratio. From Equation (23), the overall nodal coefficients \(\tilde{\mathbf{u}_{I}}\) are defined by  
\[\tilde{\mathbf{u}_{I}}(\mathbf{x})=\mathbf{u}_{h}\big{(}\mathbf{x}_{I}\big{)}-k_{1}\Bigg{[} \begin{array}{c}Q_{11}\big{(}\mathbf{x}_{I}\big{)}\\ Q_{12}\big{(}\mathbf{x}_{I}\big{)}\end{array} \Bigg{]}-k_{2}\Bigg{[} \begin{array}{ c}Q_{21}\big{(}\mathbf{x}_{I}\big{)}\\ Q_{22}\big{(}\mathbf{x}_{I}\big{)}\end{array} \Bigg{]} \tag{30}\]  
and the overall stress and strain fields corresponding to \(\tilde{\mathbf{u}_{I}}\) are recovered by the global recovery technique. Meanwhile, the original essential boundary condition (2) is modified as  
\[\hat{\mathbf{u}}=\hat{\mathbf{u}}-k_{1}\Bigg{[} \begin{array}{c}Q_{11}(\hat{\mathbf{x}}) \\ Q_{12}(\hat{\mathbf{x}})\end{array} \Bigg{]}-k_{2}\Bigg{[} \begin{array}{c}Q_{21} (\hat{\mathbf{x}})\\ Q_{22}(\hat{\mathbf{x}})\end{array} \Bigg{]},\quad\hat{\mathbf{x}}\quad on\quad\Gamma_{D} \tag{31}\]  
Substituting Equations (23) and (24) into Equation (22), together with Equations (4) and (5), ends up with  
\[\left[ \begin{array}{ccc}\mathbf{K}_{IJ}&K_{I1}&K_{I2}\\ K_{I1}^{T}&K_{I1}&K_{I2}\\ K_{I2}^{T}&K_{I2}^{T}&K_{I2}\end{array} \right]\left(