Multi Image

Triangulation from image with six connected components (thanks to F. Hecht and G. Sadaka):

Be aware that labels are not quite 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$ (this problem does not require labels)

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

 Figure Solution

load "isoline"
real hsize=5;
/* 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
*/

/*
to bluid pgm file : on shell do :
# convert filename.jpg filename.pgm
*/
real[int,int] Curves(3,1);
int[int] be(1);
int nc;// nb of curve
{
string rio;

rio="baie6.pgm";
rio="a.pgm";
rio="freefem.pgm";

real[int,int] ff1(rio); // 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 at 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; f1[]=ff1;
real vmax = f1[].max ;
real vmin = f1[].min ;
real vm = (vmin+vmax)/2;
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=vm,close=0,Curves,beginend=be,smoothing=.005,ratio=0.1);
verbosity=1;
}

int ic0=be(0), ic1=be(1)-1;

macro GG(i)
border G#i(t=0,1)
{
P=Curve(Curves,be(i*2),be(i*2+1)-1,t);
label=i+1;
}
real lg#i=Curves(2,be(i*2+1)-1);  // END OF MACRO

// For the FreeFem ++ image use these lines because there are 6 closed curves
GG(0)  GG(1)	GG(2)	GG(3)	GG(4) 	GG(5) // number of closed curves

real hh= -10;

// For the FreeFem ++ image use these lines because there are 6 closed curves
func bord = G0(lg0/hh)+G1(lg1/hh)+G2(lg2/hh)+G3(lg3/hh)+G4(lg4/hh)+G5(lg5/hh);

plot(bord,wait=0);
mesh Th=buildmesh(bord);
plot(Th,wait=0);