Let $\Omega\subset \mathbb{R}^d$ be an open bounded convex domain. Let $f\in L_2(\Omega)$. We actually know that $f$ doesn't have higher regularity. Let the pair $(\Psi,\vec{\Phi})\in L_2\times H(\mathrm{div})$ be the solution to the following elliptic PDE in mixed form.
\begin{align*} 0&= \vec{\Phi}+\vec{\nabla} \Psi& \mathrm{for\ } x\in \Omega\\ f&= \vec{\nabla}\cdot \vec{\Phi}& \mathrm{for\ } x\in \Omega\\ 0&= \Psi& \mathrm{for\ } x\in \partial\Omega_D\\ 0&= \vec{\Phi}\cdot \vec{\eta}& \mathrm{for\ } x\in \partial\Omega_N \end{align*}
We assume that the measure Dirichlet boundary $\partial\Omega_D$, in nonzero.
Does there exist a constant $C$, dependent only on the domain, such that
$$\left\| \vec{\Phi} \right\|_{L^2} + \sum_{i=1}^d \left\|\frac{\mathrm{d} }{\mathrm{d} x_i}\Phi_i \right\|_{L^2} \le C \left\| f \right\|_{L^2}?$$
My attempt
We use elliptic regularity as follows:
\begin{align*} &\left\| \vec{\Phi} \right\|_{L^2} + \sum_{i=1}^d \left\|\frac{\mathrm{d} }{\mathrm{d} x_i}\Phi_i \right\|_{L^2}\\ &\le \left\| \vec{\Phi} \right\|_{L^2} + \left| \vec{\Phi} \right|_{H^1}\\ &= \left\| \vec{\Phi} \right\|_{L^2} + \left| \Psi \right|_{H^2}\\ &= \left| \Psi \right|_{H^1} + \left| \Psi \right|_{H^2} \\ &\le \left\| \Psi \right\|_{H^2}\\ &\le C \left\|f\right\|_{L^2} \end{align*}
We use elliptic regularity in the last inequality. Note that here is where we use the fact that the domain is convex. Is there a problem with this estimate? All references I have seen for elliptic regularity results assume that there is no Neumann boundary. I appreciate any references to the literature where elliptic regularity is proven for this case.
As long as the Dirichlet boundary is of positive measure your proof looks correct. Note that the Neumann boundary condition is necessary for the well-posed part. Once you know your solution exists and is unique, you may replace Neumann with a Dirichlet condition.