Solving PDE with non-trivial null space

from netgen.geom2d import unit_square from ngsolve import * mesh = Mesh(unit_square.GenerateMesh(maxh=0.3)) h1=H1(mesh,order=5,dirichlet="bottom|left|right|top",complex=True) u,v = h1.TnT() a = BilinearForm(h1) a += SymbolicBFI(grad(u)*grad(v)) b = BilinearForm(h1) b += SymbolicBFI(u*v) a.Assemble() b.Assemble() evecs = GridFunction(h1,multidim=50) vlam = ArnoldiSolver(a.mat,b.mat,h1.FreeDofs(),evecs.vecs,1) Draw(evecs) # get a non multidim GF out of the multidim one evec0 = GridFunction(h1) evec0.vec.data = evecs.vecs[0] number = NumberSpace(mesh,complex=True) fes = FESpace([h1,number]) (u,unum), (v,vnum) = fes.TnT() # The orthogonality condition is satisfied because vnum = 1 means that evec0 * u = 0 # unum incorporates the evec0 part of the right hand side, if the argument of f is orthogonal to # evec0, then unum should be 0 alam = BilinearForm(fes) alam += SymbolicBFI(grad(u) * grad(v) - vlam[0] * u * v + unum * v * evec0) alam += SymbolicBFI(evec0 * u * vnum) f = LinearForm(fes) f += SymbolicLFI(1*v) alam.Assemble() f.Assemble() sol = GridFunction(fes,"solution") sol.vec.data = alam.mat.Inverse(fes.FreeDofs()) * f.vec Draw(sol.components[0])