load "ppm2rnm"
load "isoline"
string italy="italy.pgm";
real Areaitaly = 250000; // km^2
real [int,int] ff1(italy);
int nx = ff1.n;
int ny = ff1.m;
mesh Th1 = square(nx-1,ny-1,[(nx-1)*x,(ny-1)*(1-y)]);
fespace Vh1(Th1,P1);
Vh1 f1;
f1[] = ff1;
real [int,int] Curves(3,1);
int[int] be(1);
int nc = isoline(Th1,f1,iso=0.5,close=1,Curves,beginend=be);
int ic0 = be(0), ic1 = be(1)-1;
plot([Curves(0,ic0:ic1),Curves(1,ic0:ic1)],wait=1);
int NC = Curves(2,ic1)/10;
border G(t=0,1) 
{
  P = Curve(Curves,ic0,ic1,t);
}
plot(G(-NC),wait=1);
mesh Th = buildmesh(G(-NC));
real scale = sqrt(Areaitaly/Th.area);
Th = movemesh(Th,[x*scale,y*scale]);
cout << "Area Italy = " << Areaitaly << " km^2 " << endl;
cout << "Th.area    = " << Th.area << " km^2 " << endl;
plot(Th,wait=1);
// a priori adaptivity
fespace Vh(Th,P1);
Vh u,v;
solve Poisson(u,v)=int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
  -int2d(Th)(v)
  +on(G,u=0);
Th = adaptmesh(Th,u,iso=true,err=0.01);
plot(Th,wait=1);