load "ppm2rnm"
load "isoline"
string italy="italy2.pgm";
real Areaitaly = 250000; // km^2
real [int,int] ff1(italy); // two-dimensional image
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; // f1 is now a FE function representing Italy
real [int,int] Curves(3,1);
int[int] be(1); // actual length does not matter
int nc = isoline(Th1,f1,iso=0.5,close=1,Curves,beginend=be);
cout << "nc " << nc << endl;
int ic0 = be(0), ic1 = be(1)-1;
plot([Curves(0,ic0:ic1),Curves(1,ic0:ic1)],wait=1,cmm="curves");
int ic2 = be(2), ic3 = be(3)-1;
plot([Curves(0,ic2:ic3),Curves(1,ic2:ic3)],wait=1,cmm="curves");
int ic4 = be(4), ic5 = be(5)-1;
plot([Curves(0,ic4:ic5),Curves(1,ic4:ic5)],wait=1,cmm="curves");
int ic6 = be(6), ic7 = be(7)-1;
plot([Curves(0,ic6:ic7),Curves(1,ic6:ic7)],wait=1,cmm="curves");
// Curves(0,ic0:ic1) x coordinates first (unique in this case) curve
// Curves(1,ic0:ic1) y coordinates first (unique in this case) curve
// Curves(2,ic0:ic1) length
int NC = Curves(2,ic1)/10;
border G(t=0,1) P = Curve(Curves,ic0,ic1,t);
plot(G(-NC),wait=1,cmm="border");
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,cmm="initial mesh");
// 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.1); // Hessian of u will be used
plot(Th,wait=1,cmm="adapted mesh");
Th = adaptmesh(Th,u,iso=true,err=0.01);
plot(Th,wait=1,cmm="adapted mesh");