If you want to avoid using a Lagrange multiplier (because it enlarges the matrix to invert), some other ways to penalize the normal boundary conditions are possible.
In 3D you can use for example the boundary vertices. You have to define first a continuous approximate normal on the boundary (label 1 is for boundary)
fespace Vh1(Th,[P2,P2,P2]);
varf bdryn([w1,w2,w3],[v1,v2,v3])=int2d(Th,1)(v1*N.x+v2*N.y+v3*N.z);
Vh1 [normalappx,normalappy,normalappz];
normalappx[]=bdryn(0,Vh1);
Then you can add to your problem
+int2d(Th,qft=qf1pTlump,1)(1e10*(u1*normalappx+u2*normalappy+u3*normalappz)*(v1*normalappx+v2*normalappy+v3*normalappz))
Setting qft=qf1pTlump is a trick to get the integral on the boundary triangle performed on the three vertices of the triangle. Thus the int2d plays the role of an int0d on the boundary vertices.
Another way is to use the boundary edges:
+int2d(Th,qft=qf2pT,1)(1e10*(u1*normalappx+u2*normalappy+u3*normalappz)*(v1*normalappx+v2*normalappy+v3*normalappz))
This edge formula will work if your velocity fespace is P2, but not for P1 (there are too many edges compared to the boundary degrees of freedom of P1).
You can find the explainations in
https://perso.math.u-pem.fr/bouchut.francois/normal.pdf
Example files (2d and 3d) are:
normal2d.edp (1.3 KB)
normal3d.edp (2.1 KB)