I'll try to answer to your questions in their order of appearance below.

>1) How is the variation $\delta y$ different from the differential $dy$? one of the books i tried to use to understand this says the former is "virtual" while the later is "real". Now I am an engineer so this might sound like a dumb thing to say among mathematicians, but they both sound like they mean the same thing. maybe there is a formality here that i am not understanding. Are differentials a special case of variations ?

**Technically speaking, differentials are variations of function(al)s defined of domains on $\Bbb R$ or $\Bbb C$** (or more generally finite dimensional vector spaces), thus **they are special cases of the variation concept**. In this respect, perhaps it would also be useful to have a look at [this Q&A on functional derivatives](https://mathoverflow.net/questions/349057/question-about-functional-derivatives/349584#349584): the variation $\delta T$ of the argument of your functional is not a simple variable but a function. For example, in our case,
$$
\delta T=\varepsilon \theta \quad\varepsilon\in]0,1]
$$
where $\theta$ is a relatively small (respect to the sought for stationary $T$ value) but admissible temperature field increment, meaning with this that $I(T+\theta)$ is perfectly defined. "*Virtuality*" is an attribute given to variations that arises from the application of the [principle of virtual work](https://en.wikipedia.org/wiki/Virtual_work), itself an application of the variational principle of least action: however, they are simply variations.

>2) Why wasn't the Euler-Lagrange used here ? is it because we've now discretized the integral and it is now a summation ? or is there another reason ?

**The Euler Lagrange-Equation is not used here since it is useless in this context**: the basic aim of the method of finite elements (and likewise of the finite difference method, despite their large differences) is to approximate a vector field with with a piecewise constant vector function (a *simple function* using the terminology of integration theory) obtained by solving a linear algebraic system, hopefully sufficiently small to be tractable. In the case under analysis, you proceed by dividing your domain in a  (finite) number $n$ of regions (elements) and evaluate the given functional $I(T)$ in each of them by doing the following steps in this exact order
1. *Approximate the temperature field $T$ inside each region* by a piecewise constant function, and then
2. *Minimize the local functional*, now transformed in a simpler linear algebraic relation *by calculating its functional derivative and impose its vanishing*.

If you apply the functional derivative *before* choosing an approximation for the temperature field, you get *exactly* the Euler-Lagrange equation of the functional $I(T)$ but, as you can check, this is simply the local formulation of the former problem you are trying to solve, i.e. the stationary equation for the heat conduction
$$
\frac{\partial}{\partial x}\left(k_x\frac{\partial T}{\partial x}\right)  
+ \frac{\partial}{\partial y}\left(k_y\frac{\partial T}{\partial y}\right) + \frac{\partial}{\partial z}\left(k_z\frac{\partial T}{\partial z}\right) = G\quad \quad (x,y,z)\in \Omega_e,\; e=1, \dots, M,
$$
therefore you haven't done any step forward towards the resolution of your problem (but possibly a backward step).