Chesapeake

Triangulation from image with one connected component (thanks to F. Hecht):

Be aware that labels are not well computed and you better look at the mesh file before doing something that needs labels. Once the triangulation is obtained we solve on the image $\Omega$ (does not require labels)

\begin{array}{rcl} -\Delta u+u & =& x\,y,\quad \Omega\\ \partial_{n}u &=&0,\quad \partial \Omega. \end{array}  Triangulation Solution

real hmin=3;
/*  Chesapeake bay: Build mesh
----------------------------------
make during the European Finite Element Fair 2009
in TTK/ Helsinki  (6 june 2009)
http://math.tkk.fi/numericsyear/fefair/
By F. Hecht

FROM DATA  Chesapeake.jpg come form the  presentation of
Mme. Sonia Garcia email: smg at usna.edu

Professors Denny Kirwan and Bruce Lipphardt - UD
Mrs. Lisa Becktold and the CADIG staff
Dr. Tom Gross - NOAA
USNA Chesapeake Bay Group:
Reza Malek-Madani, Kevin McIlhany, Gary Fowler, John Pierce,
Irina Popovici, Sonia Garcia, George Nakos, Louise Wallendorf,
Bob Bruninga, Jim DArchangelo

Remark, I will all tool to build best mesh, it is a  first test
to buid mesh form image.

2 hours of work
*/

/*
to bluid pgm file : on shell do :
# convert Chesapeake.jpg Chesapeake.pgm
*/

string Chesapeake="Chesapeake.pgm";
real[int,int] ff1(Chesapeake); // 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 ll;
for(int i=0;i<pends.n;i++)
{
real xx= pends[i++];
real yy= ny-pends[i];
ll = ll + i*(square(x-xx)+(square(y-yy)) < 25);
}
plot(ll,wait=1);*/
Vh f1,f2;
f1[]=ff1; //  transforme array in finite element function.

f2= f1 < 0.99;

//plot(f1,wait=1,dim=2,fill=1);

if(hmin<6)
if(hmin<2)
//plot(Th,wait=1);

f2=f2;
Th=trunc(Th,f2>0.9); //   remove

//plot(Th,wait=1);

//  a small problem to separed the connect component of the meshes.
//  We solve a   $\epsilon u - \Delta u = 0, \partial_n u = 1$
// so the solution is almost constant on each connect component $C$
// and the value is  $1/\epsilon \int_{\partial C}/ \int_C 1$
// so the greatest component have the mininql value.

Vh u,v;
//plot(u,wait=1,fill=1);
//  set:  the  threshold value value :  l=  min + 10%
real l=u[].min + (u[].max-u[].min)*0.1;
cout << " l = " << l << endl;

Th=trunc(Th,u<l);//  remove a element greater to l
plot(Th,wait=0);

solve Problem1(u,v)=
int2d(Th)( dx(u)*dx(v)+dy(u)*dy(v)+u*v) +int2d(Th)(-x*y*v);

plot(u,cmm="Laplacian on the Image ",wait=0,fill=1);