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?