/usr/share/doc/freefem++/examples/examples++-load/Leman-mesh.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 | load "ppm2rnm" load "isoline"
string leman="lg.pgm";
real AreaLac = 580.03; // $Km^2$
real hsize= 5;
real[int,int] Curves(3,1);
int[int] be(1);
int nc;// nb of curve
{
real[int,int] ff1(leman); // read image and set to an rect. array
int nx = ff1.n, ny=ff1.m; // grey value beetween 0 to 1 (dark)
// build a Cartesian mesh such that the origne is qt the right place.
mesh Th=square(nx-1,ny-1,[(nx-1)*(x),(ny-1)*(1-y)]);
// warning the numbering is of the vertices (x,y) is
// given by $ i = x/nx + nx* y/ny $
fespace Vh(Th,P1);
Vh f1; f1[]=ff1; // transforme array in finite element function.
nc=isoline(Th,f1,iso=0.25,close=1,Curves,beginend=be,smoothing=.1,ratio=0.5);
verbosity=1;
}
// the longuest isoline
int ic0=be(0), ic1=be(1)-1;
plot([Curves(0,ic0:ic1),Curves(1,ic0:ic1)], wait=1);
int NC= Curves(2,ic1)/hsize;
real xl = Curves(0,ic0:ic1).max-5;
real yl = Curves(1,ic0:ic1).min+5;
border G(t=0,1) { P=Curve(Curves,ic0,ic1,t); label= 1 + (x>xl)*2 + (y<yl);}
plot(G(-NC),wait=1);
mesh Th=buildmesh(G(-NC));
plot(Th,wait=1);
real scale = sqrt(AreaLac/Th.area);
Th=movemesh(Th,[x*scale,y*scale]); // resize the mesh to have the correct scale
cout << " Th.area = " << Th.area << " Km^2 " << " == " << AreaLac << " Km^2 " << endl ;
plot(Th,wait=1,ps="leman.eps");
|