I have encountered an interesting problem with using Mathematica 12 to solve a nonlinear finite element problem.

Here is some sample code that generates the error:

Needs["NDSolve`FEM`"]

mesha = ToElementMesh[Rectangle[]];

\[Alpha] = 10;

F[x_, y_] = 
  Inactivate[{{D[u1[x, y], x], D[u2[x, y], x]}, {D[u1[x, y], y], 
      D[u2[x, y], y]}}] + IdentityMatrix[2];

pst = {Inactive[D][\[Alpha] (Det[F[x, y]] - 1), x], 
   Inactive[D][\[Alpha] (Det[F[x, y]] - 1), y]};

BCs = {DirichletCondition[u1[x, y] == 0, y == 0], 
   DirichletCondition[u2[x, y] == 0, y == 0], 
   DirichletCondition[u2[x, y] == .01, y == 1]};

us = NDSolveValue[{pst == {0, 0}, BCs}, {u1, u2}, {x, y} \[Element] 
    mesha];

Here, NDSolveValue generates the error: NDSolveValue::underdet: There are more dependent variables, {u1[x,y],u2[x,y],u1[x,y],u2[x,y]}, than equations, so the system is underdetermined.

What I cannot seem to figure out is why Mathematica thinks there are two unique u1's and two unique u2's. There is no typographical difference in any of the u2's or u1's typed into the code as far as I can tell. Any ideas on how to resolve this? I've tried writing the code in many different ways but continuously run into problems getting it in the proper Inactive format for Mathematica...

Thanks in advance for any advice!