macro minmax(phi,vphi) cout<<phi+" : min = "<<vphi[].min<<" ; "<<" max = "<<vphi[].max<<endl//
real dens;
dens=10;
real Lmin,Lmax,L,Hmin,Hmax;
Lmin=-1.0;Lmax=1.0;L=0.3;Hmin=2.0;Hmax=4.0;
real D,dt,usdt,v;
D=1.0;dt=0.1;usdt=1.0/dt;v=1;
border b1(s=0,1){x=Lmin+-1*L*s+-1*Lmin*s;y=0;label=1;};int nb1=max(2.0,dens*abs((2*L*Lmin+L^2+Lmin^2)^0.5));
border b2(s=0,1){x=-L;y=Hmin*s;label=1;};int nb2=max(2.0,dens*abs(Hmin));
border b3(s=0,1){x=2*L*s+-L;y=Hmin;label=1;};int nb3=max(2.0,dens*abs(2*L));
border b4(s=0,1){x=L;y=Hmin+-1*Hmin*s;label=1;};int nb4=max(2.0,dens*abs(Hmin));
border b5(s=0,1){x=L+-1*L*s+Lmax*s;y=0;label=2;};int nb5=max(2.0,dens*abs((-2*L*Lmax+L^2+Lmax^2)^0.5));
border b6(s=0,1){x=Lmax;y=Hmax*s;label=1;};int nb6=max(2.0,dens*abs(Hmax));
border b7(s=0,1){x=Lmax+-1*Lmax*s+Lmin*s;y=Hmax;label=1;};int nb7=max(2.0,dens*abs((-2*Lmax*Lmin+Lmax^2+Lmin^2)^0.5));
border b8(s=0,1){x=Lmin;y=Hmax+-1*Hmax*s;label=1;};int nb8=max(2.0,dens*abs(Hmax));
plot(b1(nb1)+b2(nb2)+b3(nb3)+b4(nb4)+b5(nb5)+b6(nb6)+b7(nb7)+b8(nb8));
mesh m=buildmesh(b1(nb1)+b2(nb2)+b3(nb3)+b4(nb4)+b5(nb5)+b6(nb6)+b7(nb7)+b8(nb8));
plot(m);
fespace mP1(m,P1);
fespace mP0(m,P0);
mP1 c1,c0,cp;
c0=0;
problem unpas(c1,cp)=int2d(m)(usdt*c1*cp+D*((dx(c1))*(dx(cp))+(dy(c1))*(dy(cp))))+int2d(m)(-1*usdt*c0*cp)+int2d(m)(1*v*(dy(cp))*c1)+on(2,c1=1);
real[int] viso=0.0:0.05:1;
real Cp;
for (int i=0;i<1000;i++) {
unpas;
c0=c1;
Cp=int2d(m)(c0);
cout << "  N = " << Cp << " ; " ;
minmax("c0",c0);
plot(c0,fill=1,viso=viso);
};