nowait;
border(1,0,1,6)  begin x:=0;      y:=1-t   end;
border(2,0,1,15) begin x:=2*t;    y:=0     end;
border(2,0,1,10) begin x:=2;      y:=-t    end;
border(2,0,1,20) begin x:=2+3*t;  y:=-1    end;
border(2,0,1,35) begin x:=5+15*t; y:=-1    end;
border(3,0,1,10) begin x:=20;     y:=-1+2*t end;
border(4,0,1,35) begin x:=5+15*(1-t); y:=1 end;
border(4,0,1,40) begin x:=5*(1-t);y:=1     end;

buildmesh(1500);

nu := 0.001; dt := 0.4;

/* initial values */
u = y*(1-y)*one(y>0);
v = 0;
p = 2*nu*x*one(y>0);
un = u; vn = v; 
i:=1; j:=2; k:=3;

iter(180)
{  
  f=convect(un,u,v,dt);  g=convect(vn,u,v,dt);
  
  solve(u,i)
    {         
      /*Horizontal velocity*/
	 onbdy(1) u = y*(1-y);
      onbdy(2,4) u = 0;
      onbdy(3)dnu(u)=0;
      pde(u) id(u)/dt-laplace(u)*nu = f/dt -dx(p);
    };
  plot(u);
  
  solve(v,j)
    {        /* Vertical velocity */
		  onbdy(1,2,3,4) v = 0;
    pde(v) id(v)/dt-laplace(v)*nu = g/dt -dy(p);
    };
  
  solve(p,k) {      /*  Pressure */
			  onbdy(1,2,4) dnu(p) = 0;
  onbdy(3) p=0;
  pde(p)  -laplace(p)= -(dx(f) + dy(g))/dt;
  };
  un = u;  vn = v; i:=-1; j:=-2; k:=-3;
}