$L^p$-norm of the gradient of a harmonic function. Let $1<p<\infty$, $\Omega \subsetneq U \subset \mathbb{R}^n$ be bounded open sets in $\mathbb{R}^n$ with smooth boundary. The claim is that there is $c>0$ such that the following inequality holds
$$ \| \nabla u \|_{L^p(\Omega^\prime)} \leq c \| u \|_{L^p(U)}$$ 
for all $u \in L^p(U)$ with $u$ being a weak solution of $\Delta u = 0$, 
I believe that it can be proved by the mean value property of harmonic functions that such a $u$ is actually smooth and so $\nabla u$ is defined in the usual sense. My question is how can we prove this inequality? The proof I'm reading just states that applying the mean value property to $\nabla u$ and so the $c$ exists. Also, I'm kind of suspicious of for the case $p \in (1,2)$ because I feel we might need to apply Jensen's inequality to $\frac{p}{2}$ somehow.
 A: In what follows, I assume that $\Omega \Subset U$. Fix $r>0$ small enough so that $B(x,r) \subseteq U$ for any $x \in \Omega$. 
Assuming that we already know the mean value property holds for $\nabla u$, for $x \in \Omega$ we have 
\begin{align*}
\nabla u(x) & = \frac{1}{|B(x,r)|} \int_{B(x,r)} \nabla u(y) dy \\ 
& = \frac{1}{|B(x,r)|} \int_{\partial B(x,r)} u(y) \vec{n}(y) dy \\ 
& = \frac{n}{r} \cdot \frac{1}{|\partial B(x,r)|} \int_{\partial B(x,r)} u(y) \vec{n}(y) dy
\end{align*}
by Gauss formula. Now we apply Jensen's inequality to obtain 
\begin{align*}
|\nabla u(x)|^p 
& \leqslant \frac{n^p}{r^p} \cdot \frac{1}{|\partial B(x,r)|} \int_{\partial B(x,r)} |u(y)|^p dy \\
& = C(n,r) \int_{\partial B(0,r)} |u(x+y)|^p dy. 
\end{align*}
Integrating over $x \in \Omega$, changing the order of integration and substituting $z = x+y$ yields 
\begin{align*}
\int_{\Omega '} |\nabla u(x)|^p dx
& \leqslant C(n,r) \int_{\partial B(0,r)} \int_{\Omega} |u(x+y)|^p dx dy \\
& \leqslant C(n,r) \int_{\partial B(0,r)} \int_{U} |u(z)|^p dz dy \\
& = C(n,r) \int_{U} |u(z)|^p dz. 
\end{align*}
