## 2 Natural Neighbour Radial Point Interpolation Method
### Nodal Connectivity and Numerical Integration
\begin{array}{ l}\mathbf{x}_{q_{l}}^{T}\\ \mathbf{x}_{q_{l}}^{T}\\ \mathbf{x}_{q_{l}}^{T}\\ \mathbf{x}_{q_{l}}^{T}\end{array} \right]=\left[ \begin{array}{l}N_{1}&N_{2} &N_{3}&N_{4}\\ N_{1}&N_{2}&N_{3}&N_{4}\\ \mathbf{x}_{q_{l}}^{T}\end{array} \right]\cdot\left[ \begin{array}{l}x_{q_{1}} &y_{q_{1}}\\ x_{q_{2}}&y_{q_{2}}\\ x_{q_{3}}&y_{q_{3}}\\ x_{q_{4}}&y_{q_{4}}\end{array} \right] \tag{2}\]  
in which \(\mathbf{N}(\xi,\eta)=\{N_{1}(\xi,\eta),N_{2}(\xi,\eta),N_{3}(\xi,\eta),N_{4}( \xi,\eta)\}^{T}=\{N_{1},N_{2},N_{3},N_{4}\}^{T}\). Being isoparametric quadrilateral with four nodes, the generic equation of \(N_{i}(\xi,\eta)\) is defined as  
\[N_{i}(\xi,\eta)=\frac{1}{4}\cdot(1+\hat{\xi}_{i}\cdot\xi)\cdot(1+\hat{\eta}_{i }\cdot\eta) \tag{3}\]  
with \(\hat{\xi}_{i}\) and \(\hat{\eta}_{i}\) being  
\[\left\{ \begin{array}{l}\Xi=\{\hat{\xi}_{1},\hat{\xi}_{2},\hat{\xi}_{3},\hat {\xi}_{4}\}=\{-1,-1,1,1\}\\ \mathbf{H}=\{\hat{\eta}_{1},\hat{\eta}_{2},\hat{\eta}_{3},\hat{\eta}_{4}\}= \{-1,1,1,-1\}\end{array} \right. \tag{4}\]  
The seminal study on the NNRPIM suggests using \(1\times 1\) or \(2\times 2\) integration points per quadrilateral sub-area [31]. The main disadvantage of using a \(2\times 2\) integration scheme (as the one represented in Figure 2) is the inherent prohibitive computational cost. In general, using a \(2\times 2\) integration scheme leads to computational analysis about \(10\sim 15\) times slower than the \(1\times 1\) integration scheme. Moreover, the gain in accuracy of a \(2\times 2\) integration scheme is not significant when compared with the \(1\times 1\) integration scheme [22,31]. If the \(1\times 1\) integration scheme per quadrilateral sub-cell is considered, the isoparametric coordinates and weight are these: \(\{\xi_{l},\eta_{l},\omega_{l}\}^{T}=\{0,0,4\}^{T}\). On the other hand, if the \(2\times 2\) integration scheme per quadrilateral sub-cell is assumed, the following isoparametric coordinates and weights must be assumed,  
\[\left\{ \begin{array}{l}\xi_{l}\\ \eta_{l}\\ \omega_{l}\end{array} \right\}=\left\{ \begin{array}{l}\frac{-1}{\sqrt{3}} &\frac{-1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}\\ \frac{-1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{-1}{\sqrt{3}} \\ 1&1&1&1\end{array} \right\} \tag{5}\]  
Subsequently, it is possible to calculate the Cartesian integration weight of each integration point with the following expression:  
\[\hat{\omega}_{l}=\frac{A_{sub}}{A_{iso}}\cdot\omega_{l} \tag{6}\]  
in which the area of the quadrilateral sub-area is defined by \(A_{sub}\) and the area of the isoparametric cell is identified as \(A_{iso}\). Notice that in this 2D case, \(A_{iso}\) is always \(A_{iso}=2\times 2=4\). Consequently, assuming a generic function \(f(x,y)\) defined inside a quadrilateral sub-area domain, \(\Omega_{s}\), it is possible to integrate \(f(x,y)\) using a Gauss-Legendre integration scheme with \(2\times 2=4\) integration points:  
\[F=\int_{\Omega_{s}}f(x,y)d\Omega_{s}=\sum_{l=1}^{4}f(x_{l},y_{l})\cdot\hat{ \omega}_{l} \tag{7}\]  
considering the Cartesian coordinates \(\mathbf{x}_{l}=\{x_{l},y_{l}\}^{T}\) and integration weights \(\hat{\omega}_{l}\) of the corresponding integration points calculated with Equations (2) and (6), respectively. This procedure is repeated for each sub-area of the Voronoi cell \(V_{i}\), and then the complete process is repeated again for each Voronoi cell of the Voronoi diagram \(\mathbf{V}\). In the end, a background set of \(n_{Q}\) integration points is obtained, allowing us to integrate any function defined within the problem's domain \(\Omega\). Since \(A_{\Omega}=\sum_{i=1}^{N}A_{V_{i}}\), thus \(A_{\Omega}=\sum_{i=1}^{n_{Q}}\hat{\omega}_{I}\). The formal construction of the Voronoi diagram and the background set of integration points is described with detail in the literature [22].