## 2 Natural Neighbour Radial Point Interpolation Method
### Discrete System of Equations  
The global system of equations for an elasticity problem can be obtained with the virtual work principle. Assuming that the work produced by external forces is equal to the work produced by the internal forces, \(W_{int}=W_{ext}\), the following expression can be written:  
\[\int_{\Omega}\delta\mathbf{\varepsilon}^{T}\cdot\mathbf{\sigma}\cdot d\Omega=\int_{ \Omega}\delta\mathbf{u}(\mathbf{x}_{I})^{T}\cdot\mathbf{b}\cdot d\Omega+\int _{\Gamma_{t}}\delta\mathbf{u}(\mathbf{x}_{I})^{T}\cdot\mathbf{\bar{t}}\cdot d\Gamma \tag{20}\]  
with the problem's domain represented by \(\Omega\). As the expression shows, on the domain's surface boundary \(\Gamma\), natural and essential boundaries, \(\Gamma_{t}\) and \(\Gamma_{u}\), respectively, can be defined. Thus, external forces \(\mathbf{\bar{t}}\) can be applied on \(\Gamma_{t}\) and displacement constrains can be imposed at \(\Gamma_{u}\). The solid domain can also be under the influence of body forces \(\mathbf{b}\), acting uniformly on \(\Omega\). Since Equation (20) requires both displacement components \(\{u,v\}\), they can be defined simultaneously:\[\mathbf{u}(\mathbf{x}_{l})= \begin{Bmatrix}u(\mathbf{x}_{l})\\ v(\mathbf{x}_{l})\end{Bmatrix} =\mathbf{H}(\mathbf{x}_{l})\cdot\mathbf{u}= \begin{bmatrix}\phi_{1}(\mathbf{x}_{l})&0&\ldots&\phi_{n}(\mathbf{x}_{l})&0\\ 0&\phi_{1}(\mathbf{x}_{l})&\ldots&0&\phi_{n}(\mathbf{x}_{l})\end{bmatrix}  \cdot\begin{Bmatrix}u_{1}\\ v_{1}\\ \vdots\\ u_{n}\\ v_{n}\end{bmatrix} \tag{21}\]  
leading to the following representation of the deformation vector:  
\[\boldsymbol{\varepsilon}(\mathbf{x}_{l})=\mathbf{L}\cdot\mathbf{u}(\mathbf{x }_{l})= \begin{bmatrix}\frac{\partial}{\partial x}&0\\ 0&\frac{\partial}{\partial y}\\ \frac{\partial}{\partial y}&\frac{\partial}{\partial x}\end{bmatrix} \cdot \Phi(\mathbf{x}_{l})\cdot\mathbf{u}=\mathbf{B}(\mathbf{x}_{l})\cdot\mathbf{u} =\mathbf{B}(\mathbf{x}_{l})\cdot\left\{ \begin{aligned} u_{1}\\ v_{1}\\ \vdots\\ u_{n}\\ v_{n}\end{aligned} \right\} \tag{22}\]  
in which the deformation matrix \(\mathbf{B}(\mathbf{x}_{l})\) is defined as  
\[\mathbf{B}(\mathbf{x}_{l})= \begin{bmatrix}\frac{\partial\phi_{1}(\mathbf{x}_ {l})}{\partial x}&0&\ldots&\frac{\partial\phi_{n}(\mathbf{x}_{l})}{ \partial x}&0\\ 0&\frac{\partial\phi_{1}(\mathbf{x}_{l})}{\partial y}&\ldots&0&\frac{ \partial\phi_{n}(\mathbf{x}_{l})}{\partial y}\\ \frac{\partial\phi_{1}(\mathbf{x}_{l})}{\partial y}&\frac{\partial\phi_{1}( \mathbf{x}_{l})}{\partial x}&\ldots&\frac{\partial\phi_{n}(\mathbf{x}_{l})}{ \partial y}&\frac{\partial\phi_{n}(\mathbf{x}_{l})}{\partial x}\end{bmatrix}  \tag{23}\]  
Considering [PERSON]'s law, it is possible to calculate the stress at the integration point \(\mathbf{x}_{l}\) with:  
\[\boldsymbol{\sigma}(\mathbf{x}_{l})=\mathbf{c}(\mathbf{x}_{l})\cdot \boldsymbol{\varepsilon}(\mathbf{x}_{l})=\mathbf{c}(\mathbf{x}_{l})\cdot \mathbf{B}(\mathbf{x}_{l})\cdot\mathbf{u} \tag{24}\]  
In a 2D problem, the material constitutive matrix \(\mathbf{c}(\mathbf{x}_{l})\) can be defined for for plane stress or plane stress conditions [47]. If plane stress is considered, the material constitutive matrix is defined with  
\[\mathbf{c}=\frac{E_{2}}{1-\alpha_{1}\cdot v_{21}^{2}} \begin{bmatrix}\alpha_{1 }&\alpha_{1}\cdot v_{21}&0\\ \alpha_{1}\cdot v_{21}&1&0\\ 0&0&\alpha_{2}\cdot(1-\alpha_{1}\cdot v_{21}^{2})\end{bmatrix}  \tag{25}\]  
If plane strain conditions are being considered, the constitutive matrix should be  
\[\mathbf{c}=\frac{E_{2}}{\alpha_{3}\cdot\alpha_{4}} \begin{bmatrix}\alpha_{1} \cdot(1-\alpha_{1}\cdot v_{21}^{2})&\alpha_{1}\cdot v_{21}\cdot\alpha_{3}&0\\ \alpha_{1}\cdot v_{21}\cdot\alpha_{3}&1-v_{12}^{2}&0\\ 0&0&\alpha_{2}\cdot\alpha_{3}\cdot\alpha_{4}\end{bmatrix}  \tag{26}\]  
in which \(\alpha_{1}=E_{1}/E_{2}\), \(\alpha_{2}=G_{12}/E_{2}\), \(\alpha_{3}=1+v_{12}\) and \(\alpha_{4}=1-v_{12}-2\cdot\alpha_{1}v_{21}^{2}\).