Well, since no one has given a proper answer I will describe the regularity for the two equations I am familiar with. Maybe someone else knows how to answer in the case of the wave equations. For the sake of simplicity, lets assume Dirichlet boundary conditions in both cases, although this does not really affect the answer.
Poisson's Equation $\Delta u = f$:
This is an elliptic PDE. The usual way to attack the regularity of these is to use the Elliptic Regularity Theorem. A very basic form of this says, for example:
If $\Omega$ is a "nice" open set in $R^n$, $f \in H^k(\Omega)$ and $u$ satisfies $\Delta u = f$, then $u \in H^{k+2}(\Omega)$.
This implies a variety of generalizations to more general equations. For example, if $\Psi \in L^{\infty}(\Omega)$, $f \in L^2(\Omega)$ and $u \in L^2(\Omega)$ satisfies $\Delta u + \Psi u = f$, then we know that $\Delta u = f - \Psi u$. But the RHS is in $L^2(\Omega)$ by the assumptions on $u$, $f$ and $\Psi$. Thus $u \in H^2(\Omega)$. Also notice that the need to know a priori that $u \in L^2(\Omega)$ is not really a restriction since it is trivial by Lax-Milgram that $u \in H^1(\Omega)$. Moreover, if $f \in C^{\infty}(\Omega)$, then $f \in H^k(\Omega)$ for every $k > 0$. Thus if $\Delta u = f$ we know that $u \in H^{k+2}$ for every $k > 0$. Then the Sobolev Embedding Theorem implies that $u$ is also smooth. In general, you should expect solutions to be classical if $f \in H^k(\Omega)$ with $k + 2$ large enough so that the Sobolev Embedding implies that $u \in C^2_b(\Omega)$ These are the type of things to keep in mind when dealing with elliptic PDE.
Heat equation $u_t = \Delta u$, $u(0, x) = f(x)$:
This is a prime example of an evolution problem (it is also parabolic). These tend to have the property that the operator on the RHS ($\Delta$ in this case) generates a continuous (sometimes even analytic) semigroup (I recommend Roger & Renardy's book for an accessible introduction). It turns out that $\Delta$ generates an analytic semigroup $e^{t\Delta}$ on $L^2(\Omega)$ so the time regularity comes pretty much for free. Moreover, if $T(t)$ is an analytic semigroup with generator $A$ then $T(t)$ maps into the domain of $A^k$ for every $k > 0$. In the case of $A = \Delta$, this means that $e^{t\Delta}$ maps into $H^k(\Omega)$ for every even $k > 0$. Thus by applying the Sobolev Embedding Theorem, $e^{t\Delta}$ maps into $C^{\infty}(\Omega)$. Therefore, for every $f\in L^2(\Omega)$, the solution $u(t) = e^{t\Delta}f$ is a strong solution on $(0, \infty) \times \Omega$ and is immediately smoothed. So for this type of equation you get both existence and regularity if you can show that the operator on the RHS generates an analytic semigroup. A similar approach works for the Schrodinger equation, which is hyperbolic, so the main thing to notice here is the "evolution form" $u_t = Au$. So for these the solution will almost always be classical
I have never studied the wave equation, so in this case I have no idea.
I can tell you my favourite proof of the mean value property, which I find more intuitive than the one via Green's identity. To make the proof rigorous, you need to know about integration over the manifold $O(n)$ of all orthogonal matrices, but you can just depict what I do as averaging over all orthogonal matrices. Let's take $x=0$ for simplicity. Heuristically, you'd expect that
$$
\frac{1}{r^{n-1}\sigma_{n-1}}\int_{S_r(0)}u(y)\,d\sigma(y) = {\rlap{\;\bar{}}\int_{O(n)}} u(Az)\,d\sigma(A) =: f(z)
$$
for $z\in\mathbb R^n$ with $|z|=r$. On the left hand side you take the average of $u(y)$ over all $y$ with $|y|=r$, on the right hand side it's the average of $u(Az)$ over all the orthogonal matrices $A$, which should be and is indeed the same.
Now you can differentiate the right hand side under the integral sign: take the Laplacian with respect to $z$. Then, due to $\Delta u=0$, you see that $f$ is harmonic. Moreover, $f$ is radially symmetric, i.e., $f(z)$ depends on $|z|$ only. Finally, use the fact that a radially symmetric harmonic function defined on all of $\mathbb R^n$ is constant. This yields $f(z) = f(0) = u(0)$, and we're done.
Best Answer
One of the cases you need it is for prescribing boundary values. Just consider the change of variable formula. To define a $k$th order derivative space on the boundary of a set, you need that boundary to be described by a function that is at least $W^{k,\infty}$ (which is the Lipschitz space $C^{k-1,1}$).
Similarly, the usual technique in proving boundary estimates for PDEs is to "straighten out" the boundary so that you can work in the case where your domain is the half-space. To do this change of variable, and preserve strong solutions in a uniform way, requires certain minimum differentiability assumptions on the boundary.
Another case where it is used is that for a bounded domain with $C^2$ boundary, the boundary is compact, and therefore as a lower bound on the radius of curvature. This is then sufficient to guarantee uniform interior and exterior sphere conditions†on the boundary, and thus allows you to prove strong maximum principle for second order elliptic operators on the domain.
†The sphere conditions are conditions on the "one sided" smoothness of the boundary. Let $\Omega$ be an open domain and let $x_0$ be a point on the boundary $\partial\Omega$. The boundary $\partial\Omega$ is said to satisfy an interior sphere condition with radius $r$ at $x_0$ if there exists $x'\in\Omega$ such that the open ball of radius $r$ around $x'$, which we denote $B_r(x')$, satisfies $B_r(x') \subset \Omega$ and $x_0 \in \partial B_r(x')$. Immediately you see that if $\partial\Omega$ satisfies the interior sphere condition with radius $r$, then it will also for all radii $r' < r$. The exterior sphere condition is defined analogously, with $x'$ and $B_r(x')$ required to reside to the exterior of $\Omega$.
Observe that the sphere conditions give "one-sided" bounds on the "radius of curvature". For example, let $\Omega = \{ (x,y)\in \mathbb{R}^2, y > |x| \}$. Then at the corner at the origin, $\Omega$ does not satisfy the interior sphere condition for any $r$, but it satisfies exterior sphere conditions for any positive radius.
Lastly, uniformity is just the statement that $\exists r_0$ st. $\forall x\in\partial\Omega$, $\partial\Omega$ satisfies an interior/exterior sphere condition of radius $r_0$ at $x$.