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