This is a problem in a PDEs course I am helping a student with, and I would like to take this opportunity to verify my solution. I will also entertain alternative/easier approaches, because I suspect I am overthinking it.
FYI: $\Delta_3=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}$
This student (and just about all his classmates) have never taken a course in analysis, so I would enjoy reading and sharing solutions that aren't too rigorous.
Here is my approach.
Suppose, for the sake of contradiction, that the max of $u$ is attained at $x_0\in D-\partial D$
Define $E=\{r>0: B_r(x_0)\subseteq D\}$
For each $r\in E$ let $\bar{u}(r)$ be the average value of $u$ on $S_r(x_0)=\{\|x-x_0\|=r\}$.
We may parameterize $S_r(x_0)$ to write $\bar{u}(r)$ in the following way: $$\bar{u}(r)=\frac{1}{4\pi r^2}\iint_{S_r(x_0)} udS=\frac{1}{4\pi r^2}\int_0^{2\pi} \int_0^\pi u(x_0+r \vec{n}(\phi,\theta))r^2 \sin(\phi)d\phi d\theta$$
Here, $\vec{n}(\phi,\theta)=(\sin(\phi)\cos(\theta),\sin(\phi)\sin(\theta),\cos(\phi))$ is your standard parameterization of the unit sphere.
Notice $\vec{n}(\phi,\theta)=\frac{x-x_0}{r}$ is an outward point normal vector to the surface $S_r(x_0)$ at $x=x_0 + r \vec{n}(\phi,\theta)$.
Also notice $u_r=\nabla u \cdot x_r$.
So use divergence theorem to get $$\begin{eqnarray*} \bar{u}'(r)&=&\frac{1}{4\pi r^2 }\int_0^{2\pi} \int_0^\pi \big(\nabla u\big)(x_0+r \vec{n}(\phi,\theta))\cdot \vec{n}(\phi,\theta)r^2 \sin(\phi)d\phi d\theta \\&=& \frac{1}{4\pi r^2} \iint_{S_r(x_0)} \Big(\nabla u\cdot \vec{n}\Big)dS \\ &=& \frac{1}{4\pi r^2} \iiint_{|x-x_0|\leq r} \text{div}\Big(\nabla u\Big)dV \\ &=& \frac{1}{4\pi r^2} \iiint_{|x-x_0|\leq r} \Delta_3 u dV \\ &>& 0\end{eqnarray*}$$ This shows $\bar{u}$ (strictly) increases with $r$.
But we also know $\bar{u}(r)\rightarrow u(x_0)$ as $r\rightarrow 0^{+}$, so let's choose $\epsilon>0$.
Find $\delta>0$ so that $\rho\in(0,\delta)$ implies $u(x_0)-\bar{u}(\rho)<\epsilon$.
Now fix $r\in (0,\delta)$.
Evidently, $$0\leq u(x_0)-\bar{u}(r/2)<\epsilon$$ $$0\leq u(x_0)-\bar{u}(r)<\epsilon$$ But this means $$\bar{u}(r)-\bar{u}(r/2)=(u(x_0)-\bar{u}(r/2))-(u(x_0)-\bar{u}(r))<\epsilon$$ We've just shown $\forall \epsilon>0: \bar{u}(r)< \epsilon + \bar{u}(r/2)$ which means $\bar{u}(r)\leq \bar{u}(r/2)$ a contradiction.