// example 4
// Navier-Stokes equations : flow around a cylinder
// 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 = 300;

border ba(t=0,1.0){x=t*5.0-1.0;y=-1.0;label=1;};
border bb(t=0,1.0){x=4.0;y=2.0*t-1.0;label=2;};
border bc(t=0,1.0){x=4.0-5.0*t;y=1.0;label=3;};
border bd(t=0,1.0){x=-1.0;y=1.0-2.0*t;label=4;};
border cc(t=0,2*pi){x=cos(t)*0.25;y=sin(t)*0.25;label=5;};
mesh Th=buildmesh(ba(30)+bb(30)+bc(30)+bd(30)+cc(-30));

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,3,u2=0)  
+ on(4,u1=1.0-y*y,u2=0) 
+ on(5,u1=0,u2=0);

real area = int2d(Th)(1.0);
real meanp;

plot(Th, wait=1);
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);

}