As far as I know, multigrid uses stationary iterative methods as smoothers (i.e GS), but can we use a direct method also?
For example, in case we have a tridiagonal system (for example 1D heat equation), can we use Thomas algorithm with spatial coarsening (in a V-cycle for example). I know that using Thomas in this specific case saves a lot of work but I am curious?
If you can solve the linear system with a direct solver, then that's exactly what you should be doing. Multigrid is a method that can be used if you don't have the time or memory resources to use a direct solver (because direct solvers have a complexity that grows faster than $O(N)$ with the size of a linear system). If you use a direct solver as a sub-step in the multigrid method, then multigrid necessarily inherits the complexity of the direct solver and will scale just as poorly.
There is one exception: In a multigrid hierarchy, at some level or other you will get to a place where the linear system becomes relatively small (because the size of the linear system is reduced by a factor of 4 (in 2d) or 8 (in 3d) every time you coarsen the mesh by one level). At that point, it becomes inefficient to continue using a multigrid method, and one just solves the linear system exactly -- typically with a direct solver. But the point is that one doesn't want to do that on the finer levels.
I don't see why you wouldn't be able to do it. So, I do not see a technical limitation.
However, using a direct method to solve (even a tridiagonal system) can be an overkill. You might be able to get away with only a few iterations (implying only a couple of very cheap matrix-vector products) to get a reasonable smoothing using an iterative method, while you would have to actually solve the system using the direct method.
Mathematically, I would express it, as follows. Define $t_\text{MVP}$ to be the time to perform a matrix-vector product, and $t_\text{Thomas}$ to solve a tridiagonal system with the same matrix.
For conventional workflow, you would require $t_\text{MVP}N_\text{iter}(\varepsilon)$ to achieve the desired accuracy $\varepsilon$. So, you would prefer to use the iterative technique when
$$ t_\text{MVP}N_\text{iter} < t_\text{Thomas} $$
which would be true in a lot of cases. Here, $N_\text{iter}(\varepsilon)$ is the number of iterations to converge to the desired accuracy. For multigrid smoothers, $\varepsilon$ is usually not too high, thus making iterative solvers a general choice.
Moreover, the more complicated is the matrix, the more costly it is to perform a direct solve.
