To rule out the trivial case that one of the $f_j$ is a nonzero constant, I interpret the linear independence condition as saying that no nontrivial linear combination of the $f_j$ is constant.

The set $\{\det J = 0\}$ is not necessarily finite. Consider for example $u_1 = r\sin(\theta)$ and $u_2 = r^2\sin(2\theta)$ in $B_1 \subset \mathbb{R}^2$, which both vanish on the $x$ axis (hence have parallel gradient there), and have linearly independent boundary data.

In dimension $n = 2$, the set $\{\det J = 0\}$ has measure zero. One can argue as follows: assume not, and let $u,\,v$ be the pair of harmonic functions. Since $u_y + iu_x$ and $v_y + iv_x$ are holomorphic, the sets $\{|\nabla u| = 0\}$ and $\{|\nabla v| = 0\}$ are locally finite, and furthermore, the imaginary part of their ratio,
$$\frac{\nabla^{\perp}u \cdot \nabla v}{|\nabla v|^2},$$
is harmonic. By hypothesis, this function vanishes on a set of positive measure, and by its harmonicity it vanishes identically in $\{|\nabla v| > 0\}$. We conclude in this set that $\nabla u = \lambda(x)\nabla v$. Harmonicity of $u$ and $v$ implies that $\nabla \lambda \cdot \nabla v = 0$, and differentiating the relation $u_x = \lambda v_x$ in $y$, the relation $u_y = \lambda v_y$ in $x$, and subtracting gives $\nabla \lambda \cdot \nabla^{\perp}v = 0$. It follows that $\lambda$ is constant in $\{|\nabla v| > 0\}$ and hence that $u - \lambda v$ is constant, contradicting linear independence of the boundary data.