Stokes.edp

//maillage

border a(t=1,0){x=0;y=0.5+0.5*t;label=1;};
border b(t=0,1){x=5*t;y=0.5;label=2;};
border c(t=1,0){x=5 ;y=0.5*t;label=2;};
border d(t=0,1){x=5+15*t;y=0;label=2;};
border e(t=0,1){x=20 ;y=t;label=3;};
border h(t=1,0){x=20*t;y=1;label=4;};
mesh Th = buildmesh( a(7) + b(40) + c(10) + d(150) + e(50) + h(100));
plot(Th,wait=1,ps="solfm.ps");
savemesh(Th,"maillage.msh");


int n=1;

fespace Xh(Th,P2);
fespace Mh(Th,P1);
Xh u2,v2;
Xh u1,v1;
Mh p,q;

int i=0;
real nu=1./100.;
real dt=0.1;
real alpha=1/dt;
real eps = 10e-8;

Xh up1,up2;
problem NS (u1,u2,p,v1,v2,q,solver=Crout,init=i) =
    int2d(Th)(
              alpha*( u1*v1 + u2*v2)
             + nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
             + dx(u2)*dx(v2) + dy(u2)*dy(v2) )
             - p*q*(0.000001)
             - p*dx(v1) - p*dy(v2)
             - dx(u1)*q - dy(u2)*q
           )
  + int2d(Th) ( -alpha*
       convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 )
  + on(1,u1=1,u2=0)
  + on(2,4,u1=0,u2=0)
  + on(3,u1=0.5,u2=0)
;

for (i=0;i<=10;i++)
 {
    up1=u1;
    up2=u2;
    NS;
     plot(coef=0.2,cmm=" [u1,u2] et p ",value=1,p,[u1,u2], ps="solfinal.ps");
 } ;
ofstream ff("uh.txt"); 
ff<<u1[];