## 2 Materials and Methods
### Variational Formulation
#### 2.2.3 Domain Discretization
\right)\frac{\partial}{\partial\eta}\left( \begin{array}[] {cc}\xi&\eta&1-\xi-\eta\end{array} \right)d\eta d\xi\] \[= \int_{0}^{1}\int_{0}^{1-\xi}\left( \begin{array}{cc}0&0&0\\ 0&1&-1\\ 0&-1&1\end{array} \right)d\eta d\xi=\left( \begin{array}{cc}0&0&0\\ 0&\frac{1}{2}&-\frac{1}{2}\\ 0&-\frac{1}{2}&\frac{1}{2}\end{array} \right);\]  
\[\left(\left(\varphi_{\widehat{i}^{\prime}}\,\frac{\partial \varphi_{\widehat{j}}}{\partial\eta}\right)_{\widehat{i}\widehat{j}}\right)_{ \widehat{i}\widehat{j}} = \int_{0}^{1}\int_{0}^{1-\xi}\left( \begin{array}{c}\xi\\ \eta&\eta\\ 1-\xi-\eta\end{array} \right)\frac{\partial}{\partial\xi}\left( \begin{array}[] {cc}\xi&\eta&1-\xi-\eta\end{array} \right)d\eta d\xi\] \[= \int_{0}^{1}\int_{0}^{1-\xi}\left( \begin{array}{cc}\xi&0&-\xi\\ \eta&0&-\eta\\ 1-\eta-\xi&0&\xi+\eta-1\end{array} \right)d\eta d\xi=\left( \begin{array}{cc} \frac{1}{6}&0&-\frac{1}{6}\\ \frac{1}{6}&0&-\frac{1}{6}\end{array} \right);\]  
\[\left(\left(\varphi_{\widehat{i}^{\prime}}\,\frac{\partial \varphi_{\widehat{j}}}{\partial\eta}\right)_{\widehat{i}\widehat{j}}\right)_{ \widehat{i}\widehat{j}} = \int_{0}^{1}\int_{0}^{1-\xi}\left( \begin{array}{c}\xi\\ \eta&\eta\\ 1-\xi-\eta\end{array} \right)\frac{\partial}{\partial\eta}\left( \begin{array}[] {cc}\xi&\eta&1-\xi-\eta\end{array} \right)d\eta d\xi\] \[= \int_{0}^{1}\int_{0}^{1-\xi}\left( \begin{array}{cc}0&\xi&-\xi\\ 0&\eta&-\eta\\ 0&1-\eta-\xi&\xi+\eta-1\end{array} \right)d\eta d\xi=\left( \begin{array}{cc}0& \frac{1}{6}&-\frac{1}{6}\\ 0&\frac{1}{6}&-\frac{1}{6}\end{array} \right).\]  
For the case of boundary elements, replacing \(\varphi_{\widehat{i}}\) and \(\varphi_{\widehat{j}}\) by (27) and (28) in the integrals to be calculated, and considering all combinations between functions, we have \[\left(\left(\varphi_{\widetilde{i}};\varphi_{\widetilde{j}}\right)_{ \widetilde{k}}\right)_{\widetilde{i}\widetilde{j}} = \int_{0}^{1}\!\!\left( \begin{array}{c}1-\xi\\ \xi\end{array} \right)\left( \begin{array}{cc}1-\xi&\xi\end{array} \right)d \xi=\int_{0}^{1}\!\left( \begin{array}{cc}(\xi-1)^{2}&-\xi(\xi-1)\\ -\xi(\xi-1)&\xi^{2}\end{array} \right)d\xi\] \[= \left( \begin{array}{cc}1&\frac{1}{3}\\ \frac{1}{6}&\frac{1}{3}\end{array} \right);\]  
i.e., similarly, for the second inner product, we have  
\[\left((1,\varphi_{\widetilde{j}})_{\widetilde{k}}\right)_{\widetilde{i} \widetilde{j}}=\int_{0}^{1}1\!\left( \begin{array}{c}1-\xi\\ \xi\end{array} \right)d\xi=\left( \begin{array}{c}\frac{1}{2}\\ \frac{1}{2}\end{array} \right)\]  
Thus, when assembling the matrices \(M\), \(W\), \(H\), and \(N\), which make up the system (21), each element corresponds to combinations of internal products between the base functions and/or their derivatives, as defined in (20). These matrices are symmetric, and the number of nodes in the finite element mesh determines their sizes. The research used the free software Gmsh to construct the mesh. Based on data made available by the Brazilian Institute of Geography and Statistics (IBGE), the domain chosen for the diffusion-advection equation was the territory limited by the state of Mato Grosso. The mesh has 5,417,224 triangles and 2,712,315 nodes, of which 7404 are border nodes and 2,704,911 are internal nodes.