## 3 Natural Element Approximation of Strain and Stress Fields
Then, the enhanced displacement approximation \(\vec{u}^{h}(x,y,z)\) in the 2D XNEM is expressed as  
\[\vec{u}^{h}(x)=\sum_{l=1}^{N}\left(\vec{a}_{l}+z\,\vec{r}_{l}\right)\,\varphi_ {l}(x,y)+\sum_{l=1}^{2}\sum_{k=1}^{4}\left(\vec{a}_{l}^{k}+z\,\vec{r}_{l}^{k} \right) \begin{cases}Q_{l}^{k}(r,\theta)\\ Q_{l}^{k}(r,\theta)\\ R_{l}^{k}(r,\theta)\end{cases}  \tag{19}\]  
where, \(\varphi_{l}(x,y)\) are the regular L/I functions [48; 49] and \(Q_{k}(r,\theta)\) and \(R_{k}(r,\theta)\) are the crack-tip singular functions defined by [45]  
\[\left\{Q_{l}^{k}(r,\theta)\right\}=\left\{\sqrt{\tau}sin\left( \begin{array}{ c}\theta\\ \theta\end{array} \right),\,\sqrt{\tau}cos\left( \begin{array}{c}\theta\\ \theta\end{array} \right),\,\sqrt{\tau}sin\left( \begin{array}{c}\theta\\ \theta\end{array} \right)sin(\theta),\,\,\sqrt{\tau}cos\left( \begin{array}{c} \theta\\ \theta\end{array} \right)sin(\theta)\right\} \tag{20}\]\[\left\{R_{I}^{k}(r,\theta)\right\}=\left\{\sqrt{r}sin\left(\widehat{\theta}\right),\ 0,\ 0,\ 0\right\} \tag{21}\]  
with \(\widehat{\theta}=\theta/2\). Where, \(\mathbf{b}_{I}=\left(\mathbf{a}_{I},\ \mathbf{r}_{I}\right)\) indicates the regular nodal displacement vector of the \(J\)-th node in the NEM grid shown in Figure 3, which consists of \(M\) Delaunay triangles and \(N\) nodes. Meanwhile, \(\mathbf{s}_{I}^{k}=\left(\mathbf{\hat{a}}_{I}^{k},\ \mathbf{\hat{r}}_{I}^{k}\right)\) denotes the singular vector of the \(k\)-th singular function at the \(I\)-th crack tip. The polar co-ordinates \(\left(r,\theta\right)_{I}(I=1,\ 2)\) are defined at two crack tips, as depicted in Figure 2.  
Substituting Equation (19) into Equations (12) and (13) leads to the enriched approximation of nonlinear in-plane strain \(e^{h}\) and transverse shear (T/S) strain \(\gamma^{h}\):  
\[ \begin{array}{l}\varepsilon^{h}=\sum\limits_{I=1}^{N}\left(\mathbf{H}_{L}+\mathbf{H }_{NL}\right)\varphi_{I}\mathbf{b}_{I}+\sum\limits_{I=1}^{2}\sum\limits_{k=1}^{4} \left(\mathbf{H}_{L}+\mathbf{H}_{NL}\right)\Phi_{I}^{k}\mathbf{s}_{I}^{k}\\ =\sum\limits_{I=1}^{N}\left(\mathbf{B}_{L}^{l}+\mathbf{B}_{NL}^{l}\right)\mathbf{b}_{I}+ \sum\limits_{I=1}^{2}\sum\limits_{k=1}^{4}\left(\mathbf{Q}_{I,L}^{k}+\mathbf{Q}_{I, NL}^{k}\right)\mathbf{s}_{I}^{k}\end{array}  \tag{22}\]  
\[\gamma^{h}=\sum\limits_{I=1}^{N}\mathbf{H}_{s}\varphi_{I}\mathbf{b}_{I}+\sum\limits_{ I=1}^{2}\sum\limits_{k=1}^{4}\mathbf{H}_{s}\Phi_{I}^{k}\mathbf{s}_{I}^{k}=\sum\limits_{ J=1}^{N}\mathbf{B}_{s}^{l}\mathbf{b}_{I}+\sum\limits_{I=1}^{2}\sum\limits_{k=1}^{4} \mathbf{Q}_{st}^{k}\mathbf{s}_{I}^{k} \tag{23}\]  
with \(\Phi_{I}^{k}=\left\{Q_{I}^{k},\ Q_{I}^{k},\ R_{I}^{k}\right\}\). Here, \(\mathbf{B}_{L}^{l}\), \(\mathbf{B}_{NL}^{l}\) and \(\mathbf{B}_{s}^{l}\) denote \(\mathbf{H}_{L}\varphi_{I}\), \(\mathbf{H}_{NL}\varphi_{I}\) and \(\mathbf{H}_{s}\varphi_{I}\), respectively similarly, \(\mathbf{Q}_{L}^{l}\), \(\mathbf{Q}_{NL}^{l}\) and \(\mathbf{Q}_{s}^{l}\) are defined by replacing \(\varphi_{I}\) with \(\Phi_{I}^{k}\). Referring to Figure 3, each L/I function is non-vanishing only over its support so that the displacement field approximated using L/I functions is continuous support-wise. In other words, the approximate displacement field belongs to \(C^{0}\)--continuous over the entire domain. Therefore, the first term on the right-hand side of Equation (23) may suffer from shear locking when the structure is thin and bending-dominated [40; 41]. This accuracy deterioration can be effectively resolved by indirectly interpolating the first term in the T/S strain (i.e., the non-singular term according to the MITC3+ shell concept [42]. The interpolation of non-singular T/S strain using this concept is summarized in Appendix A. The analytic manipulation of Equations (A1) and (A2) using Equations (13) and (16) provides the triangle-wise non-singular T/S strain given by [50]  
\[\gamma_{e}^{h}=\hat{\mathbf{B}}_{e}\mathbf{b}_{e} \tag{24}\]  
Figure 3: Uniform NEM grid and L/I functions.