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.Set((1,1)) gfu.components.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.