## 2 Natural Neighbour Radial Point Interpolation Method
### Discrete System of Equations
Notice that \(E_{i}\) is the [PERSON]'s modulus in material direction \(i\), and \(G_{ij}\) and \(v_{ij}\) are the elastic shear modulus and [PERSON]'s ratio in the material plane \(Oij\), respectively. Thus, considering how the stress and strain vectors were defined, it is possible to manipulate Equation (20) and obtain:  
\[\int_{\Omega}\delta(\mathbf{B}(\mathbf{x}_{l})\cdot\mathbf{u})^{T}\cdot( \mathbf{c}(\mathbf{x}_{l})\cdot\mathbf{B}(\mathbf{x}_{l})\cdot\mathbf{u}) \cdot d\Omega=\int_{\Omega}\delta(\mathbf{H}(\mathbf{x}_{l})\cdot\mathbf{u}) ^{T}\cdot\mathbf{b}\cdot d\Omega+\int_{\Gamma_{i}}\delta(\mathbf{H}(\mathbf{x }_{l})\cdot\mathbf{u})^{T}\cdot\mathbf{\tilde{t}}\cdot d\Gamma \tag{27}\]  
In this work, only small strains will be assumed (\(\delta\mathbf{B}(\mathbf{x}_{l})=0\) and \(\delta\mathbf{H}(\mathbf{x}_{l})=0\)), which allows us to simplify Equation (27) to:  
\[\int_{\Omega}\delta\mathbf{u}^{T}\cdot\mathbf{B}(\mathbf{x}_{l})^{T}\cdot \mathbf{c}(\mathbf{x}_{l})\cdot\mathbf{B}(\mathbf{x}_{l})\cdot\mathbf{u} \cdot d\Omega=\int_{\Omega}\delta\mathbf{u}^{T}\cdot\mathbf{H}(\mathbf{x}_{l} )^{T}\cdot\mathbf{b}\cdot d\Omega+\int_{\Gamma_{i}}\delta\mathbf{u}^{T}\cdot \mathbf{H}(\mathbf{x}_{l})^{T}\cdot\mathbf{\tilde{t}}\cdot d\Gamma \tag{28}\]  
\[\delta\mathbf{u}^{T}\int_{\Omega}\mathbf{B}(\mathbf{x}_{l})^{T}\cdot\mathbf{c}( \mathbf{x}_{l})\cdot\mathbf{B}(\mathbf{x}_{l})\cdot d\Omega\cdot\mathbf{u}= \delta\mathbf{u}^{T}\int_{\Omega}\mathbf{H}(\mathbf{x}_{l})^{T}\cdot\mathbf{b} \cdot d\Omega+\delta\mathbf{u}^{T}\int_{\Gamma_{i}}\mathbf{H}(\mathbf{x}_{l})^{T }\cdot\mathbf{\tilde{t}}\cdot d\Gamma \tag{29}\]  
allowing us to define the final system of equations of elasticity,\[\mathbf{K}\cdot\mathbf{u}=\mathbf{f}_{b}+\mathbf{f}_{t} \tag{30}\]  
and, consequently, calculate the problem's global displacement field:  
\[\mathbf{u}=\mathbf{K}^{-1}\cdot(\mathbf{f}_{b}+\mathbf{f}_{t}) \tag{31}\]  
It is possible to define the global stiffness matrix \(\mathbf{K}\) in its discretised form  
\[\mathbf{K}=\int_{\Omega}\mathbf{B}(\mathbf{x}_{I})^{T}\cdot\mathbf{c}(\mathbf{ x}_{I})\cdot\mathbf{B}(\mathbf{x}_{I})\cdot d\Omega=\sum_{I=1}^{N_{Q}}\mathbf{B}( \mathbf{x}_{I})^{T}\cdot\mathbf{c}(\mathbf{x}_{I})\cdot\mathbf{B}(\mathbf{x}_ {I})\cdot\hat{\omega}_{I}\cdot h(\mathbf{x}_{I}) \tag{32}\]  
as well as the global external force \(\mathbf{f}_{t}\) and body \(\mathbf{f}_{b}\) vectors, respectively, with:  
\[\mathbf{f}_{t}=\int_{\Gamma_{t}}\mathbf{H}(\mathbf{x}_{I})^{T}\cdot\bar{ \mathbf{t}}\cdot d\Gamma=0=\sum_{J=1}^{n_{Q}}\mathbf{H}(\mathbf{x}_{I})^{T} \cdot\bar{\mathbf{t}}\cdot\hat{\omega}_{J} \tag{33}\]  
\[\mathbf{f}_{b}=\int_{\Omega}\mathbf{H}(\mathbf{x}_{I})^{T}\cdot\mathbf{b} \cdot d\Omega=\sum_{I=1}^{N_{Q}}\mathbf{H}(\mathbf{x}_{I})^{T}\cdot\mathbf{b} \cdot\hat{\omega}_{I}\cdot h(\mathbf{x}_{I}) \tag{34}\]  
The thickness of the 2D solid at the location of integration point \(\mathbf{x}_{I}\) is defined as \(h(\mathbf{x}_{I})\). Since the NNRPIM uses interpolating shape functions, the imposition of essential and natural boundary conditions can be numerically implemented using the direct imposition method or the penalty methods.