derivation of variational forms of PDE directly from conservation form 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.
 A: Alright guess I'll have to try and do it myself then.
Assumptions: $\Omega$ is a nice enough compact domain that $C^\infty(\Omega)$ is dense in $H^1(\Omega)$, and $u$ is a solution of the conservation form of the diffusion equation -- for all smooth control volumes $\omega$,
$$-\int_{\partial\omega} k\nabla u\cdot\nu\;ds = \int_\omega f\; dx$$
where $\nu$ is the unit outward normal vector, $k$ is the (strictly positive) diffusion coefficient, and $f$ is in $L^2(\Omega)$.
We'll take the boundary conditions to be $u|_{\partial\Omega} = 0$ for simplicity; adding an extra term takes care of the Neumann or Robin cases.
Assume additionally that $u$ is in $H_0^1(\Omega)$.
Let $v$ be an arbitrary smooth function that vanishes on $\partial\Omega$.
Using the smooth coarea formula, we can write
$$\int_\Omega k\nabla u\cdot\nabla v\;dx = \int_\Omega k\nabla u\cdot \frac{\nabla v}{|\nabla v|}|\nabla v|\; dx = \int_{-\infty}^\infty\int_{\partial\{x: v(x) \ge t\}}k\nabla u\cdot\frac{\nabla v}{|\nabla v|}\;ds\;dt\ldots $$
but we also know that, if $t$ is a regular value of $v$, the unit outward normal to the hypersurface $v^{-t}(t)$ is equal to $-\nabla v/|\nabla v|$, hence
$$\ldots = -\int_{-\infty}^\infty\int_{\partial\{x: v(x) \ge t\}}k\nabla u\cdot\nu\;ds\;dt = \int_{-\infty}^\infty\int_{\{x: v(x) \ge t\}}f\;dx\;dt\ldots $$
The latter equality is because of our assumption that $u$ is a solution of the conservation form of the equations and that $v$ is smooth + Sard's theorem.
Finally the latter integral is just a rearrangement of the integral
$$\ldots = \int_\Omega f\,v\;dx$$
We can then pass to the limit for $v$ an arbitrary element of $H_0^1(\Omega)$ to say that $u$ is a solution of the variational problem
$$\int_\Omega k\nabla u\cdot\nabla v\;dx = \int_\Omega fv\; dx$$
for all $v$, QED.

I thought this would directly involve some distribution theory in order to be able to make statements like the gradient of an indicator function being a surface measure.
But the coarea formula lets you shove all that under the hood.
The same idea (+ some extra terms for initial conditions) seems like it should work for more problems as long as they can be expressed as a divergence in spacetime.
Please let me know if there are any mistakes or holes in the argument.
