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?