I am interested in showing continuity/boundedness of the weak solution to the following problem pde:

\begin{align*} 0 &= \mathbf{q} + \mathbf{\nabla}u && \quad x\in \Omega,\\ 0 &= \mathbf{\nabla} \cdot \mathbf{q} && \quad x\in \Omega,\\ 0 &= u && \quad x\in \partial \Omega_D,\\ g &= \mathbf{q}\cdot \mathbf{\eta} &&\quad x\in \partial\Omega_N. \end{align*}

How can I show that the norm to the weak solution to this problem is bounded by the Neumann Data? In other words, How can I show there exists a $C$ dependent only on the domain so that

$$ \| \mathbf{q}\|_{H^{\mathrm{div}}(\Omega)} \le C \| g \|_{H^{-1/2}(\partial\Omega_N)} $$

I would especially appreciate references to papers or books. If this is handled in any of the standard references (Grisvard, or Gilbarg and Trudinger) or the like, and I have missed it, could you tell me specifically where this is handled?

To the mods, I posted this question earlier today, on math.stackexchange. I have taken it down from there as I think it is more appropriate here.