macro regiong(m,dx1,dy1,dx2,dy2,drej,ind) int ind;{real xm=(dx1+dx2)/2,ym=(dy1+dy2)/2,x21=dx2+-1*dx1,y21=dy2+-1*dy1,d21=sqrt(x21*x21+y21*y21);ind=m(xm-drej*y21/d21,ym+drej*x21/d21).region;}//
real dens;
dens=5;
real Lmin,Lmax,L,Hmin,Hmax;
Lmin=-3.0;Lmax=3.0;L=1.0;Hmin=0.3;Hmax=1.5;
real D,dt,usdt;
D=1.0;dt=0.1;usdt=1.0/dt;
border b1(s=0,1){x=Lmin+-1*Lmin*s;y=0;label=1;};int nb1=max(2.0,dens*abs(Lmin));
border b2(s=0,1){x=Lmax*s;y=0;label=1;};int nb2=max(2.0,dens*abs(Lmax));
border b3(s=0,1){x=Lmax;y=Hmax*s;label=1;};int nb3=max(2.0,dens*abs(Hmax));
border b4(s=0,1){x=Lmax+-1*Lmax*s+L*s;y=Hmax;label=2;};int nb4=max(2.0,dens*abs((-2*L*Lmax+L^2+Lmax^2)^0.5));
border b5(s=0,1){x=L;y=Hmax+-1*Hmax*s+Hmin*s;label=1;};int nb5=max(2.0,dens*abs((-2*Hmax*Hmin+Hmax^2+Hmin^2)^0.5));
border b6(s=0,1){x=L+-1*L*s;y=Hmin;label=1;};int nb6=max(2.0,dens*abs(L));
border b7(s=0,1){x=-1*L*s;y=Hmin;label=1;};int nb7=max(2.0,dens*abs(L));
border b8(s=0,1){x=-L;y=Hmin+-1*Hmin*s+Hmax*s;label=1;};int nb8=max(2.0,dens*abs((-2*Hmax*Hmin+Hmax^2+Hmin^2)^0.5));
border b9(s=0,1){x=L*s+Lmin*s+-L;y=Hmax;label=3;};int nb9=max(2.0,dens*abs((2*L*Lmin+L^2+Lmin^2)^0.5));
border b10(s=0,1){x=Lmin;y=Hmax+-1*Hmax*s;label=1;};int nb10=max(2.0,dens*abs(Hmax));
border b11(s=0,1){x=0;y=Hmin*s;label=4;};int nb11=max(2.0,dens*abs(Hmin));
plot(b1(nb1)+b2(nb2)+b3(nb3)+b4(nb4)+b5(nb5)+b6(nb6)+b7(nb7)+b8(nb8)+b9(nb9)+b10(nb10)+b11(nb11));
mesh m=buildmesh(b1(nb1)+b2(nb2)+b3(nb3)+b4(nb4)+b5(nb5)+b6(nb6)+b7(nb7)+b8(nb8)+b9(nb9)+b10(nb10)+b11(nb11));
regiong(m,Lmax,Hmax,L,Hmax,0.5/dens,indp);
regiong(m,(L)*-1,Hmax,Lmin,Hmax,0.5/dens,inds);
plot(m);
fespace mP1(m,P1);
fespace mP0(m,P0);
mP1 c1,c0,cp;
c0=(region-indp)/(inds-indp);
plot(c0,cmm="c0");
problem unpas(c1,cp)=int2d(m)(usdt*c1*cp+D*((dx(c1))*(dx(cp))+(dy(c1))*(dy(cp))))+int2d(m)(-1*usdt*c0*cp);
real Cs,Cp;
real[int] viso=0.0:0.05:1.0;
for (int i=0;i<1000;i++) {
unpas;
c0=c1;
Cs=int2d(m,inds)(c0);
Cp=int2d(m,indp)(c0);
cout << "Ns = " << Cs << " ; Np = " << Cp << endl;
plot(c0,fill=1,viso=viso);
};