## 3 Finite Element Method
### Shape Functions  
As mentioned in Section 3.1, the region \(\Omega\) is divided into a finite number of quadrilateral elements. Two types of quadrilateral elements are used: 4-node quadrilateral elements (linear case) and 8-node quadrilateral elements (quadratic case). The main idea of the FEM is to approximate the unknown function at each discretized element using some approximation functions (also known as shape functions). Once the shape functions have been chosen, they are incorporated into the problem's formulation. There are many ways to choose a set of functions that form a vector basis on which to approximate the exact solution of the problem. Lagrange functions are the most common type of shape functions in the FEM. These functions are piecewise polynomials defined at each element \(e\). In this work, linear and quadratic Lagrange shape functions are used. Therefore, in the linear case, each element will have 4 nodes, while in the quadratic case each element will have 8 nodes. Expressions for the linear and quadratic Lagrange shape functions can be found in ref. [29]. Let \(\mathbf{N}=\{N_{1}\ N_{2}\ \cdots\ N_{n}\}\) with \(n=4,8\) be the basis of the Lagrange functions in the space where we seek to approximate the solution. The elements of the mesh that are used to discretize the region \(\Omega\) are defined by the local coordinates \((y_{1},y_{2})\). However, the shape functions are generally defined in parent coordinates \((\xi,\eta)\). Therefore, it is necessary to perform a transformation from local coordinates to parent coordinates and vice versa. The transformation from local coordinates \((y_{1},y_{2})\) to parent coordinates \((\xi,\eta)\) obeys the following mapping  
\[y_{1}(\xi,\eta) = \mathbf{N}(\xi,\eta)\mathbf{y}_{1}^{e},\] \[y_{2}(\xi,\eta) = \mathbf{N}(\xi,\eta)\mathbf{y}_{2}^{e}, \tag{18}\]  
where vectors \(\mathbf{y}_{1}^{e}=\left[y_{11},\ y_{12},\...\,\ y_{1n}\right]^{t}\) and \(\mathbf{y}_{2}^{e}=\left[y_{21},\ y_{22},\...\,\ y_{2n}\right]^{t}\) are the known nodal coordinates of the element, and \(n=4\) (\(n=8\)) has been used for the linear (quadratic) case. Using expressions (18) and the chain rule, the \(\xi\) and \(\eta\) derivatives of the \(i\)th shape function could be expressed in terms of the \(y_{1}\) and \(y_{2}\) derivatives as follows  
\[\frac{\partial N_{i}}{\partial\xi} = \frac{\partial N_{i}}{\partial y_{1}}\frac{\partial y_{1}}{ \partial\xi}+\frac{\partial N_{i}}{\partial y_{2}}\frac{\partial y_{2}}{ \partial\xi},\] \[\frac{\partial N_{i}}{\partial\eta} = \frac{\partial N_{i}}{\partial y_{1}}\frac{\partial y_{1}}{ \partial\eta}+\frac{\partial N_{i}}{\partial y_{2}}\frac{\partial y_{2}}{ \partial\eta}. \tag{19}\]  
Equation (19) can be written in a more compact form as  
\[\left[ \begin{array}{c}\frac{\partial N_{i}}{\partial\xi}\\ \frac{\partial N_{i}}{\partial\eta}\end{array} \right] = \mathbf{J}\left[ \begin{array}{c}\frac{\partial N_{i}}{\partial y _{1}}\\ \frac{\partial N_{i}}{\partial y_{2}}\end{array} \right], \tag{20}\]  
where  
\[\mathbf{J}= \begin{bmatrix}J_{11}&J_{12}\\ J_{21}&J_{22}\end{bmatrix} = \begin{bmatrix}\frac{\partial y_{1}}{\partial\xi} &\frac{\partial y_{2}}{\partial\xi}\\ \frac{\partial y_{1}}{\partial\eta}&\frac{\partial y_{2}}{\partial\eta}\end{bmatrix} , \tag{21}\]  
is the Jacobian matrix. The Jacobian matrix is shown below for an \(n\)-node element (\(n=4\) or \(n=8\))  
\[\mathbf{J}= \begin{bmatrix}\frac{\partial N_{1}}{\partial\xi}&\frac{\partial N _{2}}{\partial\xi}&...&\frac{\partial N_{n}}{\partial\xi}\\ \frac{\partial N_{1}}{\partial\eta}&\frac{\partial N_{2}}{\partial\eta}&...& \frac{\partial N_{n}}{\partial\eta}\end{bmatrix}\begin{bmatrix}y_{11}&y_{21} \\ y_{12}&y_{22}\\ \vdots&\vdots\\ y_{1n}&y_{2n}\end{bmatrix} . \tag{22}\]In general, each element of the discretization has a different Jacobian matrix, and entries of the Jacobian matrices are functions of \((\xi,\eta)\).