function [Ns,Coord,Nt,Vertex,Triangle,Index_vertex] = generate_mesh(Mx,My,Lx,Ly)

Ns	= Mx*My ;

x	= linspace(0,Lx,Mx) ;
y	= linspace(0,Ly,My) ;

[X,Y]	= meshgrid(x,y);
Coord	= [X(:) Y(:)];

Nt	= 2*(Mx-1)*(My-1);

Vertex  = reshape(1:Ns,My,Mx);

%[Vertex(:) Coord]

Triangle	= zeros(3,2,Nt);
Index_vertex	= zeros(Nt,3);

Label_tl	= Vertex(1:end-1,1:end-1);
Label_tr	= Vertex(1:end-1,2:end  );
Label_bl	= Vertex(2:end  ,1:end-1);
Label_br	= Vertex(2:end  ,2:end  );

%Let us scan the set of square, and for each square, fill to triangles entries.
for k=1:Nt/2
	Index_vertex(k,:)	= [Label_tl(k) Label_tr(k) Label_br(k)] ;
	coordinates	 	= Coord(Index_vertex(k,:),:);
	Triangle(:,:,k)  	= coordinates ;
	Index_vertex(k+Nt/2,:)	= [Label_tl(k) Label_bl(k) Label_br(k)] ;
	coordinates	 	= Coord(Index_vertex(k+Nt/2,:),:);
	Triangle(:,:,k+Nt/2)  	= coordinates ;
end


end
