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.
Does a mixed method solve this elliptic pde exactly if the source function is piecewise polynomial?
finite element;elliptic pde
The conjecture is false. Consider just the simplest case,$$  -\Delta u = 1$$in $\Omega=[0,1]^2$, and assume that $\partial u/\partial n=0$ all around (and fix the mean value of $u$ to make the problem well posed). For this case, the weak solution and the solution of the mixed formulation equals the strong solution of the Laplace equation.You can do the equivalent of a Taylor expansion around the origin (one of the corners of the domain). If the solution were polynomial, this expansion would have to terminate after finitely many terms. But you will find that it doesn't. In fact, the solution has a (rather weak) singularity at each of the corners (think something like $u \propto x^{7/4}y^{7/4}$ or similar -- a function for which a sufficiently high derivative becomes singular), something you can't get if the solution were indeed polynomial.