## 2 Natural Neighbour Radial Point Interpolation Method
### Shape Functions  
NNRPIM shape functions are constructed using the radial point interpolating technique, presented next. Assuming a 2D domain, \(\Omega\in\mathbb{R}^{2}\), in which a field function \(\mathbf{u}(\mathbf{x})\) is discretised with \(N\) nodes, \(\mathbf{X}=\{\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{N}\}\in\mathbb{ R}^{2}\), the interpolated value \(\mathbf{u}(\mathbf{x}_{I})\) of an given integration point, \(\mathbf{x}_{I}\in\mathbb{R}^{2}\wedge\mathbf{x}_{I}\
otin X\), can be calculated with:  
\[u(\mathbf{x}_{I})=\{\mathbf{r}(\mathbf{x}_{I})^{T},\mathbf{p}(\mathbf{x}_{I}) ^{T}\}\bigg{\{} \begin{matrix}\mathbf{a}\\ \mathbf{b}\end{matrix} \bigg{\}} \tag{8}\]  
in which \(\mathbf{r}(\mathbf{x}_{I})=\{r_{1}(\mathbf{x}_{I}),\ldots,r_{n}(\mathbf{x}_{I })\}^{T}\) is a radial basis function (RBF) and the polynomial basis function with \(m\) monomials is represented as \(\mathbf{p}(\mathbf{x}_{I})=\{p_{1}(\mathbf{x}_{I}),\ldots,p_{m}(\mathbf{x}_{I })\}^{T}\). Vectors \(a_{i}(\mathbf{x}_{I})\) and \(b_{j}(\mathbf{x}_{I})\) are the vectors of the non-constant coefficients of \(\mathbf{r}(\mathbf{x}_{I})\) and \(\mathbf{p}(\mathbf{x}_{I})\), respectively, defined as \(\mathbf{a}(\mathbf{x}_{I})=\{a_{1}(\mathbf{x}_{I}),\ldots,a_{n}(\mathbf{x}_{I })\}^{T}\) and \(\mathbf{b}(\mathbf{x}_{I})=\{b_{1}(\mathbf{x}_{I}),\ldots,b_{m}(\mathbf{x}_{I })\}^{T}\). Adding a polynomial basis function to the interpolation assures robustness and stability to NNRPIM shape functions. For example, including of a polynomial of order 2 confers \(C^{2}\) consistency and allows the NNRPIM to pass the standard patch test. The literature on the NNRPIM [22; 31] shows that it is enough to add a low-order polynomial basis, such as a constant polynomial basis (\(\mathbf{p}(\mathbf{x}_{I})=\{1\},m=1\)) or a linear polynomial basis (\(\mathbf{p}(\mathbf{x}_{I})=\{1,x_{I},y_{I}\}^{T},m=3\)). High-order polynomials increase the analysis's overall computational cost without improving the final NNRPIM solution. Since its early works [31], the NNRPIM uses a modified version of the initially proposed Multiquadratic Radial Basis Function (MQ-RBF) of [PERSON] [45], capable of taking into account the spatial dimension of the problem's domain:  
\[r_{j}(\mathbf{x}_{i})=r_{ij}=(d_{ij}^{2}+(\hat{\omega}_{I}\cdot c)^{2})^{p}= \left(\left(\sqrt{(x_{j}-x_{i})^{2}+(y_{j}-y_{i})^{2}}\right)^{2}+(\hat{\omega }_{I}\cdot c)^{2}\right)^{p} \tag{9}\]  
Regarding the MQ-RBF shape parameters \(c\) and \(p\), for 2D analyses, the literature recommends that \(c\) should be close to zero, but not zero, and \(p\) should be close to one, but not one [22], with \(c=10^{-4}\) and \(p=1-10^{-4}\) being generally used. Using only Equation (8) does not allow us to build a system of equations capable of providing the shape function of integration point \(\mathbf{x}_{I}\). Hence, previous works on radial point interpolators [28; 46] show that it is necessary to add an extra equation to build the required system of equations to assure a unique solution:  
\[\sum_{i=1}^{n}p_{j}(\mathbf{x}_{i})a_{i}(\mathbf{x}_{i})=0 \tag{10}\]  
in which \(j=\{1,2,\ldots,m\}\). In the end, a new equation matrix can be established joining both Equations (8) and (10):  
\[ \begin{bmatrix}\mathbf{R}&\mathbf{P}\\ \mathbf{P}^{T}&0\end{bmatrix} \bigg{\{} \begin{matrix}\mathbf{a}\\ \mathbf{b}\end{matrix} \bigg{\}}=\mathbf{G}\left\{ \begin{matrix}\mathbf{a}\\ \mathbf{b}\end{matrix} \right\}=\left\{ \begin{matrix}\mathbf{u}_{\mathrm{s}} \\ 0\end{matrix} \right\} \tag{11}\]  
The vector of the nodal parameters \(\mathbf{u}_{\mathrm{s}}\) is defined as  
\[\mathbf{u}_{\mathrm{s}}=\{u_{1},u_{2},\ldots,u_{n}\}^{T} \tag{12}\]  
and RBF matrix \(\mathbf{R}\) and polynomial basis matrix \(\mathbf{P}\) can be computed with \[\mathbf{R}= \begin{bmatrix}r_{11}&r_{12}&\ldots&r_{1n}\\ r_{21}&r_{22}&\ldots&r_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ r_{n1}&r_{n2}&\ldots&r_{nn}\end{bmatrix}  \tag{13}\]  
\[\mathbf{p}=