Solve BlockMatrix system without a preconditioner

Based on BlockMatrix Direct Solver, it seems that I am doomed for the start and will have to deal with a huge sparse matrix. However, I am having issues with this approach too.

## FE spaces
Hhsig = VectorValued(HDiv(mesh, order=0, RT=True), dim=2) # Space for σₕ
Hhu = L2(mesh, order=0, dim=2) # Space for uₕ
Hhphi = H1(mesh, order=1, dirichlet="bdry") # Space for φₕ
Hhlag = NumberSpace(mesh)

fes = Hhsig * Hhu * Hhphi * Hhlag # Product FE space

This is, however, prompting the following error

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
netgen.libngpy._meshing.NgException: Product space needs same dimension for all components

Any ideas for how to solve it?