## 2 Materials and Methods
### Variational Formulation
#### 2.2.3 Domain Discretization  
Discretizing the domain involves replacing the analysis in a continuous region with the analysis in a finite set of points. A standard method used for this is finite elements. Thus, over the domain of the equation, a mesh is built that is formed by geometric elements, which, in the case of the two-dimensional domain, can be triangles and quadrilaterals, among others. From this mesh and the appropriate choice of functions that make up the basis of the space \(\mathcal{V}^{(h)}\), it is possible to obtain spatial approximations for the solution via finite elements. The finite element technique is suitable for composing approximations to solve boundary value problems. It consists of dividing the solution domain into a finite number of simple subdomains, also called simple elements, and taking convenient functions (test functions) in each. This way, an approximate partial differential equation solution can be obtained locally for each element. All these domains and their functions are strategically brought together, thus generating the total solution [7]. The method consists of appropriately choosing the \(\varphi\) functions so that, together with their derivatives, they are non-zero in a small portion of the domain. This reduces the domain's integrals to integrals in part of the domain, where the functions and their derivatives are simultaneously different from zero, which simplifies the integration process [21]. The procedure boils down to the following:  
1. Create a mesh in \(\Omega\), defining which geometric figure is used for the finite elements;
2. Determine the degree of the approximation polynomials (which act on each element). Being appropriate to the problem, the base functions \(\mathcal{B}=\{\varphi_{1},\varphi_{2},...,\varphi_{n_{1}}\}\) can be chosen among the types piecewise linear, satisfying the following condition:\[\varphi_{i}(x_{j},y_{j})=\left\{ \begin{array}{ll}1,&\text{se }i=j\\ 0,&\text{se }i\
eq j\end{array} \right., \tag{22}\]  
with \((x_{i},y_{i})\) being the coordinates of the vertices of the geometric figures, also called mesh nodes [20]. To solve the system (21), the discretization of the region \(\Omega\) occurs using first-order triangular finite elements. The base functions are defined to satisfy the condition (22). Denote, with \(\{k_{e}\}_{e=1}^{nte}\), a finite family of \(nte\) triangles \(k_{e}\) that are two-by-two disjoint or have as intersection only one edge or a vertex, such that  
\[\Omega_{h}=\bigcup_{e=1}^{nte}k_{e}.\]  
Each mesh triangle is a simple subdomain, where the integrals are calculated locally. Set the triangle whose vertices are \((1,0)\), \((0,1)\), and \((0,0)\) as reference (\(\widehat{k}\)). This way, the functions \(\varphi\) restricted to this reference element are  
\[\varphi_{\widehat{1}}(\bar{\xi},\eta) = \bar{\xi}; \tag{23}\] \[\varphi_{\widehat{2}}(\bar{\xi},\eta) = \eta;\] (24) \[\varphi_{\widehat{3}}(\bar{\xi},\eta) = 1-\xi-\eta. \tag{25}\]  
The technique used to implement the finite element method requires the definition of an affine transformation \(T_{e}:\widehat{k}\to k_{e}\) that relates the reference element \(\widehat{k}\) to a generic element \(k_{e}\). That is, the way (23), (24) and (25) were defined, \(T_{e}\) must transform  
\[ \begin{array}{c}(1,0)\rightarrow(x_{1},y_{1});\\ (0,1)\rightarrow(x_{2},y_{2});\\ (0,0)\rightarrow(x_{3},y_{3}).\end{array} \]  
Thus, it is possible to define the following linear transformation:  
\[T_{e}(\bar{\xi},\eta)=\left( \begin{array}{cc}x_{1}-x_{3}&x_{2}-x_{3}\\ y_{1}-y_{3}&y_{2}-y_{3}\end{array} \right)\left( \begin{array}{c}\xi\\ \eta\end{array} \right)+\left( \begin{array}{c}x_{3}\\ y_{3}\end{array} \right) \tag{26}\]  
which is equivalent to  
\[ \begin{array}{c}x(\bar{\xi},\eta)=x_{3}+(x_{1}-x_{3})\bar{\xi}+(x_{2}-x_{3} )\eta\text{ and }\\ y(\bar{\xi},\eta)=y_{3}+(y_{1}-y_{3})\bar{\xi}+(y_{2}-y_{3})\eta\end{array} \]  
and whose Jacobian matrix is given by  
\[JT_{e}(\bar{\xi},\eta)=\left(