load "ppm2rnm"
load "isoline"
string leman="garda.pgm";
real AreaLac =  368; // Km^2 
real[int,int] Curves(3,1);
int[int] be(1);
int nc;
{ // build the curve file xy.txt ... 
real[int,int] ff1(leman); // read  image and set to an rect. array 
//  remark (0,0) is the upper, left corner.
int nx = ff1.n, ny=ff1.m; 
// build a cartesain mesh such that the origne is qt the right place.
mesh Th=square(nx-1,ny-1,[(nx-1)*(x),(ny-1)*(1-y)]); 
mesh Tf=Th;
// warning  the numbering is of the vertices (x,y) is 
// given by $  i = x/nx + nx* y/ny $

fespace Vh(Th,P1);

Vh f1,f2;
f1[]=ff1; //  transforme array in finite element function.
f2= min(f1, 0.25);

verbosity=3;
/*
 Usage of isoline
 the named parameter :
 iso=0.25   // value of iso 
 close=1, //   to force to have closing curve ... 
 beginend=be, // begin  and end of curve
 smoothing=.01, //   nb of smoothing process = size^ratio * 0.01 
 where  size is the size of the curve ...  
 ratio=0.5
 file="filename"
 
 ouptut:
  xx, yy  the array of point of the iso value 

a closed curve  number n is

in fortran notation the point of the curve are: 
(xx[i],yy[i], i = i0, i1) 
with :  i0=be[2*n],  i1=be[2*n+1];

*/
nc=isoline(Th,f1,iso=0.25,close=1,Curves,beginend=be,smoothing=.1,ratio=0.5); 
verbosity=1; 
}

int ic0=be(0), ic1=be(1)-1;		
plot([Curves(0,ic0:ic1),Curves(1,ic0:ic1)],wait=1);

// end smoothing the curve ....

border G(t=0,1) 
{    
  P=Curve(Curves,ic0,ic1,t);
  label=1;	
} 		
//cout << " .. "<<endl;
plot(G(-1000),wait=1); 
mesh Th=buildmesh(G(-1000));
plot(Th,wait=1);

real scale = sqrt(AreaLac/Th.area);
Th=movemesh(Th,[x*scale,y*scale]);
cout << " Th.area = " << Th.area << " Km^2 " << " == " << AreaLac <<  "   Km^2 " << endl ;

border G2(t=0,1) 
{    
  P=Curve(Curves,ic0,ic1,t);
  label=2;	
} 	
//fespace Vh(Th,P1); // set FE space
//Vh u,v; // set unknown and test function
//func f = 1;
//real error=0.1 ; // level of error
//problem Poisson(u,v,solver=CG,eps=1.0e-6) =
//int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v))
//- int2d(Th) ( f*(v) )
//+ on(1,u=200);

//for (int i=0;i<4;i++)
//{
//Poisson;
//Th=adaptmesh(Th,u,err=error,iso=1,thetamax=25); 
//error = error/2;
//} ;
//plot(u,wait=1);
fespace Vf(Th,P1);
Vf g;
g[]=[Curves(0,ic0:ic1),Curves(1,ic0:ic1)];

Th=adaptmesh(Th,g,iso=1,thetamax=25);

plot(Th,wait=1);