A lot of texts derive the variational form of a PDE as follows. First, life begins with a conservation law for the field $q$:
$$\partial_t \int_\omega G(q)\;dx + \int_{\partial\omega} F(q, \nabla q, \ldots)\cdot\nu\;ds = \int_\omega f\;dx$$
for all control volumes $\omega$, where $\nu$ is the unit outward normal vector to $\partial\omega$, $F$ is the flux function, $G$ takes $q$ to the density of some extensive quantity, and $f$ is a source term. We then apply the divergence theorem (assuming things are nice enough) to arrive at a strong form
$$\partial_t G(q) + \nabla\cdot F(q, \nabla q, \ldots) = f.$$
Then we multiply everything by a test function $v$ and push the divergence over onto $v$ using the usual tricks. So the argument goes conservation law $\rightarrow$ strong form $\rightarrow$ weak or variational form. I find the indirectness kind of annoying. Can we go directly from a conservation law to a variational form instead?
You could imagine that, for any $u$,
$$\int_\omega u\;dx = \int_\Omega u\cdot \mathbf 1_\omega\;dx$$
where $\mathbf 1_\omega$ is the indicator function of the set $\omega$. Likewise, assuming that the boundary of $\omega$ is nice enough, for some reasonable class of vector fields $F$, we can also say that
$$\int_{\partial\omega} F\cdot\nu\;ds = -\int_\Omega F\cdot\nabla\mathbf 1_\omega\;dx$$
in a distributional sense. Substituting these two relations into our original conservation law, we can then take arbitrary scalar multiples and sums to arrive at a variational form
$$\int_\Omega\left\{\partial_tG(q)\cdot v - F(q, \nabla q, \ldots)\cdot\nabla v - f\cdot v\right\}\,dx = \ldots$$
in a distributional sense for any simple function $v$, i.e. $v = \sum_i\alpha_i\mathbf{1}_{\omega_i}$ for some finite collection of control volumes $\{\omega_i\}$. (I've left off some boundary terms for brevity.) We can then take a limit as $v$ approaches a function in, say, $H^1(\Omega)$ or whatever space of test functions is most appropriate.
It's been a few years since I had measure theory or distribution theory. Can this argument be made rigorous, without an intermediate assumption that the PDE has a strong solution? I'd be perfectly happy if it worked for some common problem or class of problems; for example, the generalized Poisson equation where the conductivity coefficient has a jump discontinuity across some smooth hypersurface. I don't know that doing so has much practical use, I just find it more satisfying. Virtually every book I have on Galerkin methods or PDE theory either goes the conservation / strong / variational route or just starts with variational forms and ignores their derivation entirely.