Why is there no harmonic function on compact Riemannian manifold?

The product rule for the Laplacian of two functions is $$\triangle(fh) = f(\triangle h) + h(\triangle f) + 2\langle \nabla f,\nabla h\rangle.$$Stokes' theorem says that the integral of a divergence (hence of a Laplacian) over a compact manifold without boundary vanishes. So we may integrate the above to get $$0 = \int_M \triangle(fh)\,{\rm d}M = \int_Mf(\triangle h)\,{\rm d}M + \int_Mh(\triangle f)\,{\rm d}M + 2\int_M\langle\nabla f,\nabla h\rangle\,{\rm d}M$$If $\triangle f = 0$ and we take $h = f/2$ we obtain $$\int_M \|\nabla f\|^2\,{\rm d}M = 0 \implies \nabla f = 0 \implies f\mbox{ is constant}.$$(It goes without saying that we assume $M$ connected)

Ivo Terek's answer gives the most elegant answer when $M$ is oriented. (If un-oriented, pass to the orientation double cover.)

I'll answer directly your question about the strong maximum principle and ellipticity. Relative to a coordinate system, one has $$\Delta f=g^{ij}\frac{\partial^2f}{\partial x^i\partial x^j}-\frac{1}{2}g^{ij}g^{k\ell}\Big(\frac{\partial g_{j\ell}}{\partial x^i}+\frac{\partial g_{i\ell}}{\partial x^j}-\frac{\partial g_{ij}}{\partial x^\ell}\Big)\frac{\partial f}{\partial x^k}.$$ Treating this as a function on a bounded open subset of Euclidean space, then, uniform ellipticity would mean that there are positive constants $\lambda$ and $\Lambda$ with $\lambda I\leq[g^{ij}]\leq\Lambda I$, as $n\times n$ matrices. This is not necessarily the case, imagine for instance polar coordinates on the sphere, so that $g=d\theta^2+\sin^2\theta\,d\varphi^2$ and so, in this case, the upper bound $\Lambda$ does not exist. One could also imagine examples where a positive lower bound $\lambda$ does not exist.

However, the restriction of this function to any precompact subregion of the coordinate chart ensures the existence (and finiteness and positivity) of the constants, by compactness and continuity of the matrix-valued function $[g^{ij}]$.

So all you need to accept is the following statement: there is a coordinate cover $\{(U_\alpha,\varphi_\alpha)\}$ of the manifold and for each $\alpha$ an open precompact set $K_\alpha\subset U_\alpha$ such that $\{K_\alpha\}$ remains an open cover of the manifold. Relative to any coordinate chart $(K_\alpha,\varphi_\alpha)$, the operator $f\mapsto \Delta f$ is uniformly elliptic, and so one can replay the proof of the Hopf maximum principle and recover the usual strong maximum principle.

This leads to a result slightly stronger than in Ivo Terek's answer, since there the manifold must be compact for the integration to work. Here one gets the strong maximum principle, where it only has to be assumed that a maximum (or minimum) is achieved.

Also of interest is S.-T. Yau's Liouville theorem (Indiana Univ. Math. J. 1976) which says that if $(M,g)$ is a complete Riemannian manifold and $u\geq 0$ is harmonic with $\int_M u^p\,d\mu_g<\infty$ for some $p>2$, then $du=0$.