Proving the mean value property for harmonic functions using Green's second identity Suppose that $\phi:\mathbb R^3\to\mathbb R$ is a harmonic function. I am asked to show that for any sphere centered at the origin, the average value of $\phi$ over the sphere is equal to $\phi(0)$. I am also directed to use Green's second identity: for any smooth functions $f,g:\mathbb R^3\to \mathbb R$, and any sphere $S$ enclosing a volume $V$, $$\int_S (f\nabla g-g\nabla f)\cdot dS=\int_V (f\nabla^2g-g\nabla^2 f)dV.$$
Here is what I have tried: left $f=\phi$ and $g(\mathbf r)=|\mathbf r|$ (distance from the origin).  Then $\nabla g=\mathbf{\hat r}$, $\nabla^2 g=\frac1r$, and $\nabla^2 f=0$. Note also that $\int_Sg\nabla f\cdot dS=r\int_S\nabla f\cdot dS=0.$ Then $$\frac{1}{4\pi r^2}\int_S f(x)\,dA=\frac{1}{4\pi r^2}\int_S f(x)\nabla g(x)\cdot dS=\frac{1}{4\pi r^2}\int_S (f\nabla g-g\nabla f)\cdot dS.$$ Using Green's identity, this is equal to $$\frac{1}{4\pi r^2}\int_V (f\nabla^2g-g\nabla^2 f)\,dV=\frac{1}{4\pi r^2}\int_V\frac{f}{r}\,dV.$$ This reminds me of the Cauchy integral formula. If there is indeed some sort of identity that I can use to show that the last integral is equal to $f(0)$? Or is there another way to solve this problem?

I appreciate the other answers, but I came up with my own answer which, in my humble opinion, is a bit simpler.
For posterity's sake, here is the answer: Let $g=1/|\mathbf r|$ and $f=\phi$. Then the following hold: $$\int_S g(\mathbf r)\nabla f(\mathbf r)\cdot d\mathbf S=-\frac{1}{r}\int_S \nabla f\cdot d\mathbf S=0,$$ $$\nabla^2 f(\mathbf r)=0\quad\quad\text{for } \mathbf r\in V,$$ $$\nabla g(\mathbf r)=-\frac{\mathbf{\hat r}}{r^2},$$ $$\nabla^2 g(\mathbf r)=-4\pi \delta(\mathbf r).$$ Thus, the average value of $\phi$ on $S$ is
    \begin{align*}
    \frac{1}{4\pi R^2}\int_S \phi(\mathbf r)\,dS&=\frac{1}{4\pi}\int_S f(\mathbf r)\left(\frac{\hat r}{R^2}\right)\cdot d\mathbf S\\
    &=-\frac{1}{4\pi}\int_S f(\mathbf r)\nabla g(\mathbf r)\cdot d\mathbf S\\
    &=-\frac{1}{4\pi}\int_S (f(\mathbf r)\nabla g(\mathbf r)-g(\mathbf r)\nabla f(\mathbf r)\cdot d\mathbf S\\
    &=-\frac{1}{4\pi}\int_V(f(\mathbf r)\nabla^2 g(\mathbf r)-g(\mathbf r)\nabla^2 f(\mathbf r))\,dV\\
    &=\frac{4\pi}{4\pi}\int_V\phi(\mathbf r) \delta(\mathbf r)\,dV\\
    &=\phi(\mathbf 0).
\end{align*}
