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.
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*}