Using Galerkin method for PDE with Neumann boundary condition? I am wanting to show existence of solutions to
$$u_t +L(u) = f \;\;\text{on}\;\; \Omega$$
with initial condition $u|_{t=0} = u_0$ and Neumann boundary condition $\nabla u\cdot \nu = 0$ on ${\partial\Omega}$. 
How do I do this with the Galerkin method? My problem is with the boundary condition. Recall that the Galerkin method requires a triple $V \subset H$ with continuous dense inclusion. In this case, 
$$V = \{u \in H^1 : \nabla u \cdot \nu = 0 \;\;\text{on}\;\; \partial\Omega\}$$
and $H = L^2$. But showing $V$ is dense in $H$ is a problem. 
So I guess this is not how we do it. How else can I tackle this problem without putting the BC in the Hilbert space?
Thanks
 A: Let us assume the equation is second order parabolic and $L = -\Delta$ for a moment. Multiplying an arbitrary test function $v\in H^1(\Omega)$, the integration by parts of $\displaystyle -\int_{\Omega} \Delta u v\, dx$ reads:
$$
-\int_{\Omega} \Delta u v\, dx = \int_{\Omega} \nabla u \cdot \nabla v\, dx - \int_{\partial \Omega} (\nabla u \cdot \nu) v\, dS
$$
Notice that the integration on boundary does not appear in the Galerkin weak formulation:
$$
\int_{\Omega} u_t v\, dx + \int_{\Omega} \nabla u \cdot \nabla v\, dx = \int_{\Omega} f v\, dx
\tag{1}
$$
i.e., for any test function the boundary integral simply vanishes. For we didn't impose any boundary conditions on the test function $v$, $\nabla u \cdot \nu$ must vanish a.e. on the boundary. If it were not the case, we could always find a $v$ such that (1) does not hold. 

EDIT: Now suppose we solved (1) for all $v\in H^1$, the question is why $\nabla u \cdot \nu$ vanishes on boundary. If we assume $u\in H^2$ and perform integration by parts on the $\displaystyle\int_{\Omega} \nabla u \cdot \nabla v\, dx$ backwards we have
$$
\int_{\Omega} u_t v\, dx - \int_{\Omega} \Delta u \, v\, dx + \int_{\partial \Omega} (\nabla u \cdot \nu) v\, dS= \int_{\Omega} f v\, dx
\tag{2}
$$
Since $v$ is an arbitrary test function in $H^1$ (here is the trick, we just manipulating $v$ to get the results we want), we could let $v=0$ on $\partial\Omega$, then the boundary integral vanishes, and
$$
\int_{\Omega} (u_t- \Delta u)v\, dx= \int_{\Omega} f v\, dx
$$
holds for any $v\in H^1_0$. We could argue in $\Omega$: $u_t- \Delta u = f$. Now back to (2) we have $\displaystyle \int_{\partial \Omega} (\nabla u \cdot \nu) v\, dS = 0$ for any $v\in H^1$ when $v$ is not zero on boundary. Hence $\nabla u \cdot \nu = 0$.
