Finite Element Discretization of Differential Operator (FEM Programming)

Question

I am pretty new to FEM. I want to try to discretize a PDE operator using FEM---in order to find approximations of the eigenvalues and eigenfunctions for example. I know I can try to use NDSolve or NDEigensystem, but I want to check what I do by hand using mathematica (in order to try to get a deeper understanding of how the numerical methods do and how to implement them). I was hoping someone could at least point me in the right direction.

Consider a PDE operator $\mathcal L$, for example take some advection operator $\mathcal L:\phi \mapsto f^\intercal\nabla\phi$, where $f$ is some vector field $f:\in\mathbb{R}^2\to\mathbb{R}^2$. Now, I want to compute a discretization of the operator $\mathcal L$ over some domain $D=[a,b]\times[a,b]$ with $(0,0)\in[a,b]\times[a,b]$. I want to take the simplest approximation, that is, a piecewise linear approximation over a triangular mesh. To construct the triangular mesh we divide the domain $D$ into $N^2$ squares where each square has a side length $h = \frac{(b-a)}{N}$ and we further divide each square into two triangles with sides $h$ and hypotenuse $\sqrt{2}h$. So, the mesh has $2N^2$ triangular elements $T_{k}$. My aim is to programmatically construct the matrix $$ L_{i,j} = \int_{(x,y)\in T_{i}} f(x,y)^\intercal\nabla\phi_{i}(x,y)\phi_{j}(x,y){\rm d}x{\rm d}y $$ with $i,j\in\{1,\ldots,2N^2\}$. Here, $\phi_{k}$ denotes the global linear basis function on the $k$-th triangle $T_{k}$. From here, I believe I can take the left eigenvectors $l$ of $L$ to construct approximations of the eigenfunctions $g$ of $\mathcal L$, via $g(x,y) = \sum_{i=1}^{2N^2}l_{i}\phi_{i}(x,y)$.

Is there any way of doing this using the lower-level FEM functionality of Mathematica 12?. Any advice, or direction would be appreciated.