// example 3
// Navier-Stokes equations : regularized cavity flow problem 
// COE-tutorail 2007, Atsushi Suzuki 2007/03/08
//

int i;
real nu = 1.0/400.0;
real dt = 0.1;
real alpha = 1.0/dt;
int timestepmax = 200;

mesh Th=square(20,20);
fespace Vh(Th,P2),Qh(Th,P1);

Vh u1,u2,v1,v2,up1,up2;
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;

problem NS([u1,u2,p], [v1,v2,q],solver=Crout,init=i) = 
  int2d(Th)( 
           alpha * (u1*v1 + u2*v2)
         + 2.0*nu * (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
  ) 
- int2d(Th)(
           alpha * convect([up1,up2],-dt,up1)*v1
         + alpha * convect([up1,up2],-dt,up2)*v2
)
+  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;

u1 = 0;
u2 = 0;
for ( i = 0; i < timestepmax; i++) { 
   up1=u1;
   up2=u2;
   NS;
   meanp = int2d(Th)(p) / area;
   p = p - meanp;
   if (i%10 == 0) 
     plot([u1,u2],p,wait=1,value=true,coef=0.1);

}