Greetings to the NGSolve community,
I am currently using ngsolve for some numerical experiments in my PhD thesis. So far, I am quite fond of its possibilities and the nice python interface. Thanks to the programmers for this!
I constantly try to enlargen my knowledge about FEM and numerics in general, but am still at its very beginnings, so I hope this question is not too stupid.
Debugging my code I boiled some problems I had down to the following snippet:
from ngsolve import *
from netgen.geom2d import SplineGeometry
geo = SplineGeometry()
geo.AddRectangle( (0, 0), (2, 0.41), bcs = ("wall", "outlet", "wall", "inlet"))
mesh = Mesh( geo.GenerateMesh(maxh=0.05))
V = VectorH1(mesh, order=3, dirichlet="wall|outlet|inlet")
Q = H1(mesh, order=2)
X = FESpace([V,Q])
u,p = X.TrialFunction()
v,q = X.TestFunction()
a = BilinearForm(X)
a += SymbolicBFI(InnerProduct(grad(u),grad(v))+div(u)*q+div(v)*p)
a.Assemble()
gfu = GridFunction(X)
gfu.components[0].Set((1,1))
gfu.components[1].Set(1)
fdofs = X.FreeDofs()
for i in range(len(gfu.vec)):
if fdofs[i] == False:
gfu.vec[i] = 0
print("gfu.vec ", min(gfu.vec), max(gfu.vec))
t1 = gfu.vec.CreateVector()
t2 = GridFunction(X)
t1.data = a.mat * gfu.vec
t2.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="umfpack") * t1
print("t2 ", min(t2.vec), max(t2.vec))
As you can see, I am considering the Stokes problem with no-slip boundary conditions everywhere
and check how good the stiffness matrix is inverted. As a rather basic
test I compare the minimum and maximum value of t2.vec and gfu.vec and expect them to be the same up to some (small) deviations. This works perfectly if I leave out the inlet or outlet boundary in the declaration of the Dirichlet boundary. Considering full no-slip boundary, however, the inverse is done quite poorly.
I am confused by this result and up to my knowledge I am not considering an ill-posed problem,
so in theory this should work. Therefore, it seems like I am using ngsolve wrong, but cannot figure out what is missing.
I have to admin that I am not using the latest version of ngsolve. Mine is
NGSolve-6.2.1902-78-g34df28fd if this clarifies anything.
I will be quite grateful for any hint.
Best regards,
Philipp