Let $g$ be a smooth Riemannian metric on the closed $n$-dimensional unit disk $\mathbb D^n$.
Let $f: \mathbb D^n \to \mathbb{R}^n$ be a smooth orientation-preserving immersion, and let $\omega :\mathbb D^n \to \mathbb{R}^n$ be the harmonic map corresponding to the Dirichlet problem imposed by $f$, i.e. $\omega|_{\partial \mathbb D^n}=f|_{\partial \mathbb D^n}$.
Must $d\omega$ be invertible outside a set of measure zero in $\mathbb D^n$?
Here is a complete solution for the case where $g$ is real-analytic:
First, we note that $d\omega$ is singular on as set of positive measure if and only if it is singular everywhere on $\mathbb D^n$: Since $\omega$ must be real-analytic, so is $\det d\omega$, and the zero-set of a real-analytic function which is not identically zero has measure zero.
However $d\omega$ cannot be singular everywhere, since
$$ \int_{\mathbb D^n} \det d\omega = \int_{\mathbb D^n} \det df>0,$$
which finishes the job.
I am not sure how to deal with the case where $g$ is smooth but not real-analytic:
One idea I had is to approximate the smooth $g$ with analytic metrics $g(t)$, and to look at the corresponding harmonic maps $\omega(t)$: If $g(t)$ converges to $g$ in a sufficiently strong sense, then the $\omega(t)$ should converge to the harmonic map $\omega$ associated with the original metric $g$. However, I am not sure that this gets us anywhere, since the limiting map may have many more singular points compared to the $\omega(t)$. (The rank can fall in the limit).
Other observations:
We know that near $\partial \mathbb{D}^n$ $\text{rank}(d\omega) \ge n-1$ because $f$ is an immersion. Perhaps there is hope that somewhere near the boundary, $d\omega$ will be invertible.
What happens if we assume that $\text{rank}(d\omega) = n-1$ on a neighbourhood of $\partial \mathbb D^n$? Does this lead to a contradiction?
The following example shows that, at least on all $\partial \mathbb D^n$, the rank can be exactly $n-1$, if we only consider the "local problem": Set
$$\omega(r,\theta)=((r+r^{-1})\cos \theta,(r+r^{-1})\sin \theta)).$$ $\omega$ is harmonic on $\mathbb D^2 \setminus \{(0,0)\}$ and agrees with the immersion $(2x_1,2x_2)$ on $\partial \mathbb D^2$, but its rank on $\partial \mathbb D^2$ is $1$ everywhere.
$\omega$ is not really a counter-example to the question even on $\mathbb D^2 \setminus \{(0,0)\}$, since $d\omega$ is invertible in the interior.