a:= 6; b:= 1; c:=0.5; 
/*changewait;*/
border(1,0,1,8)  begin x:= -a;       y:= 1+b - 2*(1+b)*t  end;
border(2,0,1,26) begin x:= -a+2*a*t; y:= -1-b*(x/a)*(x/a)*(3-2*abs(x)/a ) end;
border(3,0,1,8)  begin x:= a;   y:=-1-b + (1+ b )* t   end;
border(4,0,1,20)  begin x:= a - a*t;   y:=0  end;
border(4,0,pi,8)  begin x:= -c*sin(t)/2; y:=c/2-c*cos(t)/2  end;
border(4,0,1,30)  begin x:= a*t;   y:=c end;
border(3,0,1,8)  begin x:= a;   y:=c + (1+ b-c )* t end;
border(5,0,1,35) begin x:= a-2*a*t; y:= 1+b*(x/a)*(x/a)*(3-2*abs(x)/a) end;
buildmesh(800);

solve(psi){
onbdy(1) psi = y;
onbdy(2) psi =-1 -b;
onbdy(3) dnu(psi) =0;
onbdy(4) psi = 0;
onbdy(5) psi =1+b;
pde(psi) -laplace(psi) =0;
};
plot(psi);

u1 = dy(psi); u2 = -dx(psi); dt =0.2; nu=0.01;
solve(w,1){
onbdy(1,2,4,5)w=0;
onbdy(3)dnu(w)=0;
pde(w) id(w)*2/dt - laplace(w)*nu = exp(-10*((x+a*0.8)^2 + (y^2)));
};
phi =w;
plot(phi);
iter(20) 
begin
  f =  convect(phi,u1,u2,dt)/dt;
  plot(f);
solve(w,-1){
onbdy(1,2,4,5)w=0;
onbdy(3)dnu(w)=0;
pde(w) id(w) * 2/dt -laplace(w) * nu 
      =  phi/dt +f;
};
  phi = 2 * w - phi;
  plot(phi);
end;