It is well known that solutions to the Laplace equation in a region $\Omega\subseteq\mathbb{R}^n$, $\nabla^2u=0$ satisfy the mean value property, namely for all $x\in\Omega$, and for all $r>0$,
$$ u(x) = \frac{1}{\left|B(x,r)\right|}\int\limits_{B(x,r)}u\ dV = \frac{1}{\left|\partial B(x,r)\right|}\int\limits_{\partial B(x,r)}u\ dS $$
where $B(x,r)$ is a ball of radius $r$ around $x$ (we assume $B(x,r)\subseteq \Omega$).
I was wondering whether there is a generalized property for solutions to the Poisson equation, namely suppose $\nabla^2u=f$ in some region $\Omega\subseteq\mathbb{R}^n$.
Can we find some function $w$ in $\Omega$ such that for all $x\in\Omega$ and for all $r>0$, the following is true?
$$ u(x) = \frac{\int\limits_{B(x,r)}uw\ dV}{\int\limits_{B(x,r)}w\ dV} = \frac{\int\limits_{\partial B(x,r)}uw\ dS}{\int\limits_{\partial B(x,r)}w\ dS} $$
The only thing I know is that we are constrained to the fact that if $\nabla^2u \geq 0$, then
$$ u(x) \leq \frac{1}{\left|B(x,r)\right|}\int\limits_{B(x,r)}u\ dV = \frac{1}{\left|\partial B(x,r)\right|}\int\limits_{\partial B(x,r)}u\ dS $$ and the other way around if $\nabla^2u\leq 0$.
Note that in your wish equality would imply that if $u$ is constant on a sphere, it is also constant on the whole ball. This is certainly not true.
I think what you're looking for is the Feynman-Kac formula, which says $$ u(x)=\frac{1}{|\partial B(x,r)|}\int_{\partial B(x,r)}u \;dS+\int_{B(x,r)}f\;d\mu, $$ where the measure $\mu$ can be expressed by a Brownian motion $W$ and its exit time $\tau$ from the ball of radius $r$: $$ \int_{B(x,r)}f\;d\mu=\mathbb{E}\int_0^{\tau}f(x+W_t)\;dt. $$