// example 1
// Stokes equations : regularized cavity flow problem 
// COE-tutorail 2007, Atsushi Suzuki 2007/03/08
//
mesh Th=square(20,20);
fespace Vh(Th,P2),Qh(Th,P1);
//fespace Vh(Th,P1b),Qh(Th,P1);

Vh u1,u2,v1,v2;
Qh p,q;
macro d11(u1)     dx(u1) //
macro d22(u2)     dy(u2) //
macro d12(u1,u2)  (dy(u1) + dx(u2))/2.0 //

real epsln = 1.0e-6;

solve stokes([u1,u2,p], [v1,v2,q]) = 
int2d(Th)( 2.0*(d11(u1)*d11(v1)+2.0*d12(u1,u2)*d12(v1,v2)+d22(u2)*d22(v2))
	    - dx(v1)*p - dy(v2)*p - dx(u1)*q - dy(u2)*q
    	    - p*q*epsln
) 
+  on(1,2,4,u1=0,u2=0)  
+  on(3,u1=x*(1-x)*4,u2=0) ;

real area = int2d(Th)(1.0);
real meanp = int2d(Th)(p) / area;
cout << "mean pressure = " << meanp << endl;
p = p - meanp;

plot([u1,u2],p,wait=1,value=true,coef=0.1);