1bool fast=true;
2load "mat_dervieux";  // external module in C++ must be loaded
3border a(t=0, 2*pi){ x = cos(t); y = sin(t);  }
4mesh th = buildmesh(a(100));
5fespace Vh(th,P1);
6
7Vh vh,vold,u1 = y, u2 = -x;
8Vh v = exp(-10*((x-0.3)^2 +(y-0.3)^2)), vWall=0, rhs =0;
9
10real dt = 0.025;
11// qf1pTlump means mass lumping is used
12problem  FVM(v,vh) = int2d(th,qft=qf1pTlump)(v*vh/dt)
13                  - int2d(th,qft=qf1pTlump)(vold*vh/dt)
14      + int1d(th,a)(((u1*N.x+u2*N.y)<0)*(u1*N.x+u2*N.y)*vWall*vh)
15+ rhs[] ;
16
17
18matrix A;
19MatUpWind1(A,th,vold,[u1,u2]);
20if(fast)
21  {
22    varf  vFVM(v,vh) = int2d(th,qft=qf1pTlump)(v*vh/dt)
23      - int1d(th,a)(((u1*N.x+u2*N.y)<0)*(u1*N.x+u2*N.y)*vWall*vh)      ;
24    real[int] rhs0=vFVM(0,Vh);
25    matrix M=vFVM(Vh,Vh,solver=CG);
26    A=-A+M;
27
28    for ( real t=0; t< pi ; t+=dt)
29      {
30	vold=v;
31        rhs[]=rhs0;
32	rhs[] += A * vold[] ;
33	v[]= M^-1*rhs[];
34	plot(v,wait=0);
35      }
36  }
37else
38for ( real t=0; t< pi ; t+=dt){
39  vold=v;
40  rhs[] = A * vold[] ;
41  FVM;
42  plot(v,wait=0);
43};
44