Question: I am reading Temam's book Navier Stokes Equations and he defines $E(\Omega) = \left\{u \in L^2\left(\Omega\right), \ \operatorname{div}(u) \in L^2\left(\Omega\right)\right\}$. Later he says that if $p \in H^1(\Omega)$ is the weak solution to the Neumann problem 
$$\Delta p = 0 \\
\frac{\partial p}{\partial n} = \phi
$$ 
and if $u = \operatorname{grad} p$ then $u \in E(\Omega)$.

I don't understand the definition of $E(\Omega)$. I thought it meant that the $i^\text{th}$ component of $u$ has an $i^\text{th}$ weak partial derivative and the sum is $0$. But how can $u$ have any more weak derivatives if it is the gradient of $p \in H^1$?

Answer: 
Let $u$ be a vector field in $L^2(\Omega)^n$. Then a function $q\in L^1\_{loc}(\Omega)$ is the divergence (in the distributional sense) of $u$ if it satisfies 
$$
\int\_\Omega u\cdot \nabla \phi = -\int\_\Omega q \phi
$$
for all functions $\phi\in C\_0^\infty(\Omega)$. This is written $q = \textrm{div} u$. Thus, the definition is similar to the definition of a weak derivative. Note, that we do not require that the individual derivatives $\frac\partial{\partial x\_i}u\_i$ exist.

Then the space $E(\Omega)$ is the space of $L^2$ vector fields whose divergence exists (in the distributional sense) and belongs to $L^2$. 

Let now $p$ be the weak solution of $-\Delta p=0$ with Neumann boundary conditions. Set $u=\nabla p$. Then $u\in E(\Omega)$: To see this take $\phi\in C\_0^\infty(\Omega)$
$$
\int\_\Omega u \cdot\nabla \phi = \int\_\Omega \nabla p \cdot \nabla \phi = 0,
$$
which implies $\textrm{div}u=0$. This motivates that divergence-free fields are gradient fields.
I'm just going to put this up here for my own reference. If someone thinks it's not correct please explain. Thanks.

I think he means div $u$ as a derivative in the distributional sense. That is we define $\left\langle \operatorname{div} \nabla p, \phi\right\rangle := -\left \langle\nabla p,\nabla \phi\right\rangle$ for all smooth compactly supported $\phi$. Furthermore, we define $\|T\|\_2 = \sup\left\{\left\langle T,\phi\right\rangle : \ \|\phi\|\_{L^2} \leq 1\right\}$ for any distribution $T$. 

Since $p$ is the solution of the Neumann equation above we have that $\int \nabla p \nabla v = \int \phi v$ for all $v \in H^1(\Omega)$ where the second integral is over the boundary of the domain. Finally, since $H^1$ includes the test functions and the right side of the last equation is $0$ for all test functions we have that $\|\operatorname{div} u\|\_2 = 0$ and thus $\operatorname{div} u = 0$ and $ \operatorname{div} u \in E$.