Intuition behind the Maximum Principle for Harmonic Functions In the text "Theory of Functions of One Complex Variables", I'm having trouble gaining intuition behind the Maximum Principle, as follows:

If $u:U \rightarrow \mathbb R$ is harmonic on a connected open set $U$ and if there is a point $P_{o} \in U$ with the property that $u(P_{o})= \sup_{Q \in U}u(Q)$, then $u$ is constant on $U$.

 A: Starting from the top: let's work in a bounded connected open set $U$ of $\mathbb{R}^n$ for simplicity. A harmonic function $u : U \to \mathbb{R}$ is a twice continuously differentiable function that satisfies Laplace's equation $\Delta u = 0$. Since $$\Delta u := \nabla \cdot \nabla u$$ this means the divergence of the gradient of $u$ is $0$. Physically speaking, when we take the sum of the rates of change of $u$ in the directions specified by, say, the standard basis vectors, we get $0$. 
Now, fix a point $x \in U$ and $r< \textrm{dist}(U,x)$. Then by the discussion above, the value $u(x)$ should be the average of the values of $u$ in a ball of radius $r$ around $x$. This is known as the $\textbf{mean value property}$: $$u(x) = \frac{1}{\textrm{vol}(B(x,r))} \int_{ B(x,r)} u(y) dy$$ 
Onto the maximum principle: suppose there is a point $a$ in the interior of $U$ at which $u$ achieves its maximum. If $u$ is not constant on $U$ there is a point $b$ with $u(b) < u(a)$, one may use continuity to find $$u(a) = \frac{1}{\textrm{vol}(B(a,r))} \int_{ B(a,r)} u(y) dy < u(a)$$ (for suitable $r>0$); contradiction. To gain more intuition for this argument, imagine you have a finite set of numbers $x_1,\dots, x_n$. If $$\max_j x_j = \frac{1}{n}\sum_{k=1}^n x_k$$ and we have a $k$ such that $x_k < \max x_j$, then $\max x_j$ can by no means be a genuine maximum, which is nonsense!
A: Here's two pieces of intuition, one mathematical and one physical.
Physical
The equation that describes the flow of heat is famously
$$ \partial_t u = \Delta u, $$
where $u=u(t,x)$ and $\Delta$ is the Laplacian with respect to $x$. Therefore Laplace's equation describes a steady situation, where heat does not flow. So suppose now that $u$ is harmonic and nonconstant, i.e. $\Delta u=0$, and that $u$ has a local maximum at $q$. Then there is a region around $q$ where $u$ has values smaller than $u(q)$. But we expect heat to flow from hotter to colder place, which suggests that heat should flow out of $q$ into the surrounding region, contradicting that $u$ is a steady state.
Mathematical
The above may look a bit suspicious, so let's look at an alternative explanation using the mean value formula. Let $u$ be harmonic, and let $D(a,r)$ be a disk centred at $a$ with radius $r$. Then the Divergence Theorem says
$$ 0 = \int_{R<|x-a|<r} \Delta u \, dV = \int_{\partial D(a,r)} \partial_r u \, dS = \frac{d}{dr} \int_{\partial D(a,r)} u \, dS $$
We see that this is therefore independent of $r$. Dividing by $2\pi r$ and taking $r \to 0$ gives the mean value formula,
$$ u(a) = \frac{1}{2\pi r}\int_{\partial D(a,r)} u(x) \, dS. $$
So if $a$ is a local maximum, $u(a)> u(x)$ for all $x$ in $\partial D(a,r)$ for $D(a,r)$, which is a contradiction since the mean integral must lie between the supremum and infimum on the circle, which is strictly smaller than $u(a)$.
