The n-dimensional Poisson problem is to solve $\Delta _n u = - f$ in a region $R_n $ with $u = g$ on $\partial R_n $. For a given ${\bf x}_ * $ in $R_n $ approximations \[ u\left( {{\bf x}_ * } \right) \simeq \sum _{j = 1}^M {A_j g\left( {{\bf x}_j } \right) + \sum _{k = 1}^N {B_k f\left( {{\bf \xi} _k } \right)} } \] are discussed. Given that the ${\bf x}_j $, $A_j $ are a harmonic interpolation formula of degree d for the Dirichlet problem for $R_n $ one wishes to find the ${\bf \xi} _k $, $B_k $ so that this approximation is exact for all polynomials of degree $ \leqq d$. Existence of the ${\bf \xi} _k $, $B_k $ is related to the question of whether or not $E[u] \equiv u({\bf x}_ * ) - \sum\nolimits_{j = 1}^M {A_j g({\bf x}_j )} $ is a positive linear functional of f of degree $d - 2$. This is established in certain cases. Some specific formulas are given as examples.