Why weak formulations in numerical mathematics? Regard the Poisson equation on the domain $\Omega = [-1,1]^n$ with $f \in H^{-1}$
$- \triangle u = f$
with homogenous Neumann boundary conditions. From standard regularity theory we know $u \in H^1$. Let us now take a look on a formulation of the problem in numerical analysis. For any test function $v \in H^1$ we have
$-\int_{\Omega} \operatorname{div} \operatorname{grad} u \cdot v = -\int_{\Omega} \operatorname{div} (\operatorname{grad} u \cdot v) + \int_{\Omega} \operatorname{grad} u \cdot (\operatorname{grad} v)$.
$= -\int_{\partial\Omega} n \circ (\operatorname{grad} u \cdot v) + \int_{\Omega} \operatorname{grad} u \cdot (\operatorname{grad} v) = \int_{\Omega} \operatorname{grad} u \cdot (\operatorname{grad} v)$
In the common Galerkin method we choose a series of finite dimensional subspaces that are defined based on a triangultation of the the domain $\Omega$. These may be the nodewise hat functions. We then look for a linear combination of hat functions, such that the above equation holds.
But we do we use the weak formulation anyways? Having shrinked the degrees of freedom to assert that our discrete solution $u_h$ fulfills the boundary conditions, there does not seem to persist a reason why to focus on the weak formulation anymore. But so is the case in literature.
A problem be that a direct discretization of
$Lu = f$
including boundary conditions might be an over-determined system. Why is this direct approach not treated (even if only to be discarded) in literature?

(I hope I did not totally misunderstand the question.) There are methods  (point collocation methods) based on the strong formulation satisfying $Lu=f$ on a finite set of points, without the need for meshing, hence called meshfree (or meshless); you find many papers and books on them. G.R. Liu's An Introduction to Meshfree Methods and Their Programming for instance describes everything in a general framwork of "weighted residuals", of which FEM, meshless methods and others are special cases (depending on the weighting function and on how are the residuals defined).
My personal experience was that many were marketing meshfree methods because of their alleged simplicity, but depending on the nature of problem you solve, additional problems may show up eventually; in my case, there were difficulties prescribing boundary conditions, or the necessity of locally determining spatial point configuration (i.e. local mesh).
