mesh Th=square(40,40,flags=3);

fespace Vh(Th,P1nc);
fespace Pdh(Th,P0);
fespace Pch(Th,P1);
fespace Ph(Th,P1dc);
macro grad(u) [dx(u),dy(u)]//
macro Grad(u) [grad(u#1),grad(u#2)]//
macro div(u) (dx(u#1)+dy(u#2))//

Vh u1,u2,v1,v2;
Pdh pd,qd;
Pch pc,qc;
real eps = 1e-8; 
solve Stokes( [u1,u2,pc,pd],[v1,v2,qc,qd])
= int2d(Th) ( (Grad(u):Grad(v)) - pd*div(v) - qd*div(u) + [v1,v2]'*grad(pc) + [u1,u2]'*grad(qc)
   -eps*pc*qc -eps*pd*qd)
+  on(3,u1=(1-x)*x,u2=0)   
+ on(1,2,4,u1=0,u2=0);

plot([u1,u2],wait=1);

Ph p = pd+pc;
plot(p, wait=1);
