Leman-mesh.edp

load "ppm2rnm" 
load "isolineP1" 
bool wait=false; 
string leman="lg.pgm"; 
real AreaLac =  580.03; // Km^2 
int nbsmooth = 200; 
int hsize= 1; 

{ // 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); 
fespace Vf(Th,P1); 
Vh f1,f2; 
Vh fu,fv; 
f1[]=ff1; //  transforme array in finite element function. 

f2= min(f1, 0.25); 

real[int] v=isolineP1(Th,"xy.txt",f1,iso=0.25); 
} 

int bcurve,ecurve; 
{ // Analyse the file to get the longest curve .... 
  int i=0,p=0,il=0,k=0,i0=0,i1=0; 
  ifstream f("xy.txt"); 
  string l; 
  while ( f.good()  ) 
    { 
      getline(f,l); 
      if( l=="") { 
	il = i-p; 
	if(il > i1-i0) {i1=i;i0=p;} ; 
	p=i  ; 
      } 
       else  ++i; 
    } 
  if(il > i1-i0) {i1=i+1;i0=p;} ; 
  cout << "The  longest  curve  = " << i0 << " " << i1 << " " << i1-i0 << endl; 
  bcurve=i0; 
  ecurve=i1; 
} 
real[int,int] Curve(2,ecurve-bcurve); 
{ 
  ifstream f("xy.txt"); 
  real px,py; 
  int k=0,i=0; 
  while (i<=ecurve) 
    { 
      i++; 
      real px,py; 
      f >> px >> py; 
      //cout << i << px << " " << py << endl; 
      if(  bcurve < i   && i <= ecurve) 
	{ 
	  Curve(0,k)=px; 
	  Curve(1,k)=py; 
	 k++; 
	} 

    } 
} 

//  Smoothing the curve 
{ 
  real[int,int] CurveO(Curve); 
  int nC=Curve.m; 
  // lissage 
  real[int,int] P(Curve),PP(Curve); 
  real c1=1, c0=4, ct=2*c1+c0; 
  c1/=ct; 
  c0/=ct; 
  for(int step=0;step<nbsmooth ;++step) 
    { 
      for(int j=0;j<nC;++j) 
      { 
	int j0= (j+ nC-1)%nC; 
	int j1= (j+ 1)%nC; 
	P(0,j) = c0*PP(0,j) +c1*PP(0,j0) +c1*PP(0,j1); 
	P(1,j) = c0*PP(1,j) +c1*PP(1,j0) +c1*PP(1,j1); 
      } 
      PP=P; 
    } 
  Curve=PP; 
  plot([Curve(0,:),Curve(1,:)],[CurveO(0,:),CurveO(1,:)], wait=1); 
} 
// end smoothing the curve .... 

border G(t=0,Curve.m) 
{ 
  int n= Curve.m; 
  int i=int(t), i1= (i+1)%n; 
  real t1=t-i, t0=1-t1; 
  i = i%n; 
  assert(t1>=0 && t1<1.); 
  assert(i>=0 && i < n ) ; 
  assert(i1>=0 && i1 < n ) ; 

  x = Curve(0,i)*t0 + Curve(0,i1)*t1; 
  y = Curve(1,i)*t0 + Curve(1,i1)*t1; 
  label=1; 
} 

plot(G(100),wait=1); 
mesh Th=buildmesh(G(500)); 
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 ; 
int Thnbvx=AreaLac/(hsize*hsize/5.)+1000; 

Th=adaptmesh(Th,hsize,IsMetric=1,nbvx=Thnbvx,errg=5); 
Th=adaptmesh(Th,hsize,IsMetric=1,nbvx=Thnbvx,errg=5); 
Th=adaptmesh(Th,hsize,IsMetric=1,nbvx=Thnbvx,errg=5); 
Th=adaptmesh(Th,hsize,IsMetric=1,nbvx=Thnbvx,errg=5); 
 // end of build 2d mesh 

plot(Th,wait=1);