Stokes equation with purely Neumann boundary conditions

I tried to implement the Stokes equation with purely Neumann boundary conditions. Since only Neumann conditions are applied, u is only determined up to a constant c. Then I added a Lagrange multiplier to the weak form so that the variational problem is well-posed.

[code]
k = 2
V = VectorH1(mesh,order=k)
Q = H1(mesh,order=k-1)
R = NumberSpace(mesh)
X = FESpace([V, Q, R, R])
gfu = GridFunction(X)

(u,p, lam1, lam2), (v,q, mu1, mu2) = X.TnT()

lam = CoefficientFunction((lam1, lam2))
mu = CoefficientFunction((mu1, mu2))[/code]

When running the code, the order of convergence of p was good, but the that of u was 0. The error of u is still the same when refining mesh. Please see the attached file.

Could you please tell me how to fix this error? Any help would be appreciated.
Thank you so much.

Attachment: Untitled9.ipynb

Hi dong,

by using the NumberSpaces you force that your solution has zero mean value for both velocity components.
Your exact solution, however, does not have zero mean value.

If you change your exact solution to have zero mean

u1=sin(pi*x)*sin(pi*y)-4/pi**2

you should get convergence (see attached file)

https://ngsolve.org/media/kunena/attachments/889/stokes_neumann.py

Best
Michael

Attachment: stokes_neumann.py

Thank you so much for your clarification.