I want to prove that a solution of the 2d problem
Find $$B \in H_0(div)=\{B \in [L^2(\Omega)]^2 | div B \in L^2(\Omega), B\cdot n=0 \ \ on\ \ \partial \Omega\}\\
E \in H^1_0(\Omega)\ \ s.t\ -\nabla div B + curl E=0$$
satisfies $\|div B\|_{L^2}=0$.
The weak form is to find $C \in H_0(div)$ s.t.
$(div B, div C)+ (curl E, C)=0$. (1)
On convex domains we can consider the equation $(div \nabla \varphi=)\Delta \varphi = div B$. Since $div B$ is in $L²$ and $\Omega$ is convex we know that $\varphi \in H^2$. Then we can test (1) with $C=\nabla \varphi$ and get $(div B, div B) + (curl E, \nabla \varphi)=0$. The second term is zero since we can do integration by parts and $curl \nabla \varphi=0$. Here we have used that $\nabla \varphi$ has enough regularity to build the curl of it.
But what do we do on a non-convex domain? Then the result that $\varphi\in H^2(\Omega)$ doesn't hold anymore. One could say $\varphi \in H^1(\Omega)$ and $\nabla \varphi \in H(div)$. But then we don't have enough regularity to perform integration by parts in the second term.
Let $\phi$ be the weak solution of $$ \Delta \phi = div \ B $$ subject to Neumann boundary conditions. That is $$ \int \nabla \phi \cdot \nabla u =- \int div \ B \cdot u = \int B \cdot \nabla u $$ for all $u\in H^1(\Omega)$. This proves that $div(\nabla \phi)=B$ in the weak sense and $\nabla \phi \cdot n =0$ in the weak sense. So $\nabla\phi \in H_0(div)$.