// -*- c++ -*-
vector n = (20,20,4);
vertex a = (-1,-1,-1);
vertex b = ( 1, 1, 1);

mesh M = structured(n,a,b);
M = periodic(M,0:1,2:3,4:5);

function u = 1;
function v = 0;
function w = 0;

femfunction phi(M) = one(x>=0.4)*one(x<=0.8)*one(y>=-0.2)*one(y<=0.2);

double i=0;
double N = 20; // number of iterations
mesh T = tetrahedrize(M);
save(medit, "phi",phi,T);
save(medit, "phi",T);

double dt = 3/N;
do {
   solve(psi) in M
    {
      pde(psi)
	-div(0.1*dt*grad(psi))+psi=convect(u,v,w, dt, phi);
     }
   
  i++;
  save(medit,"phi".i,psi,T);
  phi=psi;
} while (i<N);