Let $\Omega\subset \mathbb{R}^d$, $d\in \{2,3\}$ be an open bounded polygonal/polyhedral set. Suppose I want to solve the following pde
\begin{align*} \vec{q}+\vec{\nabla}u &=0\,&x&\in \Omega\\ \vec{\nabla} \cdot \vec{q} &= f\, \quad& x&\in \Omega\\ u&=0\,& x&\in \partial\Omega_{\mathrm{dir}}\\ \vec{q}\cdot \vec{\eta}&= 0\,& x&\in \partial\Omega_{\mathrm{neu}} \end{align*}
To avoid complications assume that the measure of the Dirichlet boundary is nonzero.
Suppose further that there is triangulation of $\Omega$, $T$, such that the source function $f$ is piecewise polynomial with respect to $T$. Suppose that $T$ has no "hanging nodes". Suppose I want to approximate the solution to the pde using a the same triangulation $T$.
If I choose polynomial trial/test spaces sufficiently rich, can I generate the exact solution to the pde? For example, if $f$ is only piecewise constant, can I choose the first Raviart Thomas space that contains (continuous) piecewise affine functions for the trail/test space for $\vec{q}$? and an appropriate trail/test space for $u$ and generate the exact solution? This result is true for $d=1$, but I haven't been able to prove it to myself for $d=2$ or higher.
Since this seems true and I haven't been able to prove it, I am looking for a counter example, a proof, or at least a sketch of a proof. I appreciate any ideas. I am mainly interested in proofs/ideas that allow for domains that are not convex.