//example31.edp from file chap3/stokes.edp
verbose = 0;
int n = 3;
mesh Th=square(10*n, 10*n);

fespace Uh(Th, P1b); 
Uh u,v,uu,vv;
fespace Ph(Th, P1);  
Ph p,pp;

solve stokes([u, v, p],[uu, vv, pp]) =
    int2d(Th)(dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv)
            + dx(p)*uu + dy(p)*vv + pp*(dx(u)+dy(v))
            -1e-10*p*pp)            
            + on(1,2,4, u=0, v=0) + on(3, u=1, v=0);

plot([u,v],p,wait=true);