I'm writing this mostly for my own benefit, since I'm currently working on the same things.
First, we show that if $h(z)$ is harmonic in a simply connected domain $D$ (a domain being a connected open set, simply connected meaning 'no holes'), then there is $F(z)$ analytic (synonymous with holomorphic; use 'real analytic' for a real-valued funcion) with $\Re F(z) = h(z)$.
As the previous poster pointed out, the first observation is that for harmonic $h(z)$, the function $\psi = \frac{\partial h}{\partial x} - i\frac{\partial h}{\partial y}$ is also harmonic. Again, as mentioned, this follows from the C.R. equations, but one must also note that we have continuity of the partials (there are many examples where the C.R. equations hold, yet the function is not analytic).
Now, since $D$ is simply connected and the function $\psi$ is analytic, it has a primitive in $D$ (simple connectivity is crucial here, as this is where Cauchy's theorem holds). Call this primitive $G(z)$, so
$$G'(z) = \frac{\partial h}{\partial x} - i\frac{\partial h}{\partial y}.\tag{1}$$
But we know that $G$ is analytic, so $G = u+iv$. Thus,
\begin{gather*}
G' = u_x + iv_x = v_y-iu_y
\end{gather*}
Comparing this last equation to (1) (using $u_x+iv_x$ for the first comparison and $v_y-iu_y$ for the second), we find
\begin{gather*}
\Re \frac{\partial G}{\partial x} = \frac{\partial h}{\partial x} = \frac{\partial \Re G}{\partial x} \cr
\Re \frac{\partial G}{\partial y} = \frac{\partial h}{\partial y} = \frac{\partial \Re G}{\partial y}
\end{gather*}
Thus,
\begin{gather*}
\frac{\partial}{\partial x}(\Re G(z) - h(z)) = 0 \cr
\frac{\partial}{\partial y}(\Re G(z) - h(z)) = 0
\end{gather*}
Since $D$ is a connected domain, we must have $h(z) = \Re G(z) + c$, so setting $F(z) = G(z) + c$ gives us $h(z) = \Re F(z)$, as desired.
We now have, given a harmonic function $h$, an analytic function whose real part is $h$. This required some work and wasn't totally trivial; in particular, we needed simple connectivity of the domain. What's important for us is the total converse: given $h(z) = \Re F(z)$ for $F$ analytic, $h$ is harmonic. You can prove this using the Cauchy-Riemann equations
Now, we want to show that the composition of a harmonic function with an analytic one is harmonic. Suppose then $\phi(z)$ is harmonic in $V$, and let $f(z)$ be analytic in $U$ with $f(U) \subset V$. Then $u(f(z))$ is harmonic for $z \in U$.
To prove this, keep in mind that harmonicity is a local property. By the previous result, we have (locally) some function $\varphi(z)$, analytic, such that $\phi(z) = \Re \varphi(z)$. Let then $z_0 \in U$ with $f(z_0) = w_0 \in V$. Then $\varphi(f(z)) = \varphi(w)$ will be analytic for $w$ near $w_0$, so
$$u(f(z)) = \Re \varphi(f(z)). $$
By the total converse, $u(f(z))$ is harmonic for $z$ close enough to $z_0$.
The following inequality helps: if $D$ is a disk centered at $z_0$, $f$ is holomorphic in $D$ and continuous on $\overline D$, then
$$\sup_{D'} |f|\le |\operatorname{Im}f(z_0)| + 3\max_{\partial D}|\operatorname{Re}f|, \tag1$$
where $D'$ is the disk concentric with $D$ and half the radius.
Proof of (1) is immediate from the Schwarz integral formula:
$$f(z)= i\operatorname{Im}f(z_0)+\frac{1}{2\pi i}\int_{\partial D}\frac{\zeta+z-2z_0}{\zeta-z} \,\operatorname{Re}f(\zeta)\,\frac{d\zeta}{\zeta-z_0}, \tag2$$
where $$\left|\frac{\zeta+z-2z_0}{\zeta-z}\right|\le 3,\qquad \zeta\in\partial D,\quad z\in D'. \tag3$$
Returning to the problem at hand, an application of (1) to the differences $f_n-f_m$ yields uniform convergence of $f_n$ on any disk $D$ centered at $z_0$ and compactly contained in $\Omega$. Now (1) can be applied to disks whose centers lie in the set in which uniform convergence is already known. Repeating in this fashion, we can cover any compact subset of $\Omega$ in finitely many steps.
Best Answer
This answer does not use complex analysis per se, but a basic property of the Lebesgue integral, the dominated convergence theorem. Sorry! I thought you might find it useful nonetheless.
To show that $u$ is harmonic, choose a point $z_0\in U$ and observe that each $u_n$ satisfies the mean-value property at $z_0$: $$ u_n(z_0) = \frac{1}{2\pi r}\int_{|z-z_0|=r} u_n(z), \text{ for all sufficiently small $r>0$.} $$ Since each $u_n$ is dominated in modulus by $1$, the dominated convergence theorem allows us to let $n$ approach infinity and interchange limit with integral. Then $$ u(z_0) = \frac{1}{2\pi r}\int_{|z-z_0|=r} u(z), \text{ for all sufficiently small $r>0$.} $$