(This really should be a comment, but is a bit long.)

Consider elliptic regularity and maximum principle. Interior regularity of harmonic function means that away from the boundary, your function $u_f$ is real analytic, and hence its level sets will also be analytic arcs away from where $\nabla u_f = 0$. Near the boundary the behaviour may degenerate, and the regularity will depend on the regularity of $f$.

Where $\nabla u_f = 0$, by the maximum principle for harmonic functions, the Hessian $\nabla^2 u_f$, if non-zero, must be indefinite (or you can see that just by noting it is trace free). So the non-degenerate critical points of $u_f$ are saddle points, and have the classical structure with two level-lines intersecting there.

In general, near a critical point, taking the Taylor expansion of the function $u_f$, you must have

$$ u_f(x) = u_f(x_0) + \sum_{|\alpha| \geq m} a_\alpha (x-x_0)^\alpha $$

where $\alpha$ are multi-indices, $a_\alpha$ are coefficients, and $m \geq 1$ is the highest number of derivatives to which $u_f$ vanish. The maximum principle states that

$$ \sum_{|\alpha| = m+1} a_\alpha x^\alpha $$

cannot be a signed function. So the critical point there correspond to an intersection of at least 2 and up to $m+1$ level curves.