bool fast=true;
load "mat_dervieux";  // external module in C++ must be loaded
border a(t=0, 2*pi){ x = cos(t); y = sin(t);  }
mesh th = buildmesh(a(100));
fespace Vh(th,P1);

Vh vh,vold,u1 = y, u2 = -x;
Vh v = exp(-10*((x-0.3)^2 +(y-0.3)^2)), vWall=0, rhs =0;

real dt = 0.025;
// qf1pTlump means mass lumping is used
problem  FVM(v,vh) = int2d(th,qft=qf1pTlump)(v*vh/dt)
                  - int2d(th,qft=qf1pTlump)(vold*vh/dt)
      + int1d(th,a)(((u1*N.x+u2*N.y)<0)*(u1*N.x+u2*N.y)*vWall*vh)
+ rhs[] ;


matrix A;
MatUpWind1(A,th,vold,[u1,u2]);
if(fast)
  {
    varf  vFVM(v,vh) = int2d(th,qft=qf1pTlump)(v*vh/dt)
      - int1d(th,a)(((u1*N.x+u2*N.y)<0)*(u1*N.x+u2*N.y)*vWall*vh)      ;
    real[int] rhs0=vFVM(0,Vh);
    matrix M=vFVM(Vh,Vh,solver=CG);
    A=-A+M;

    for ( real t=0; t< pi ; t+=dt)
      {
	vold=v;
        rhs[]=rhs0;
	rhs[] += A * vold[] ;
	v[]= M^-1*rhs[];
	plot(v,wait=0);
      }
  }
else
for ( real t=0; t< pi ; t+=dt){
  vold=v;
  rhs[] = A * vold[] ;
  FVM;
  plot(v,wait=0);
};
