I read through the Documentation and the NGS build-in solvers but I didn't figure out how to use/adapt things.
Code:
from netgen.geom2d import unit_square
from ngsolve import *
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
V = H1(mesh, order= 1)
u,v=V.TnT()
a=BilinearForm(V,symmetric=True)
a+=grad(u)*grad(v)*dx
a.Assemble()
f=LinearForm(V)
f.Assemble()
f+=v*dx
N=NumberSpace(mesh)
p,q=N.TnT()
b=BilinearForm(trialspace=V,testspace=N)
b+=u*q*dx
b.Assemble()
g=LinearForm(N)
g.Assemble()
K=BlockMatrix( [[a.mat,b.mat.T],[b.mat,None]] )
rhs=BlockVector( [f.vec,g.vec] )
gfu=GridFunction(V)
gfN=GridFunction(N)
sol=BlockVector( [gfu.vec,gfN.vec] )
sol.data=K.Inverse(inverse="pardiso")*rhs #???
