function [nodes,elem,edges]=CreateRectMesh(x,y,blow,bup,bleft,bright)
#[...] = CreateRectMesh(...)
# generate a rectangular mesh
#
# [nodes,elem,edges] = CreateRectMesh(x,y,blow,bup,bleft,bright)
# nodes elem edges  contain the information on the mesh
#
# x y  vectors containing the coordinates of the nodes
#   a typical vertex is (x(j), y(i))
#
# blow,bup,bleft,bright indicate the type of boundary condition
#     at lower, upper, right and left edge of rectangle
#   b*=1  Dirichlet condition
#   b*=2  Neumann condition
n=length(x);
m=length(y);

if (nargin!=6) usage("[nodes,elem,edges]=CreateRectMesh(x,y,blow,bup,bleft,bright)"); endif

%===================== nodes =======================================
nodes=zeros(n*m,3);
for in=1:n
  for im=1:m
    nodes((im-1)*n+in,:)=[x(in),y(im),0];
  endfor 
endfor

for im=1:m
  nodes((im-1)*n+1,3)=bleft;
  nodes(  im*n    ,3)=bright;
endfor  

if (bleft==1)  instart = 2;   else instart=1; endif 
if (bright==1) inend=(n-1);   else inend=n;   endif

for in=instart:inend
  nodes(in,    3)=blow;
  nodes((m-1)*n+in,3)=bup;
endfor

%======================= elements =====================================

elem=zeros((n-1)*(m-1)*2,4);

cc=1;
for row=1:m-1
  for col=1:n-1
    elem(cc,:)=[(row-1)*n+col,row*n+col+1,row*n+col,1];
    cc++;
    elem(cc,:)=[(row-1)*n+col,(row-1)*n+col+1,row*n+col+1,1];
    cc++;
  endfor
endfor

%========================= edges ===================================

edges=zeros(2*(n+m-2),3);
cc=1;
for in=1:n-1
  edges(cc,1:3)=[in,in+1,blow];
  cc++;
endfor

for in=1:n-1
  edges(cc,1:3)=[(m-1)*n+in,(m-1)*n+in+1,bup];
  cc++;
endfor

for im=1:m-1
  edges(cc,1:3)=[(im-1)*n+1,im*n+1,bleft];
  cc++;
endfor

for im=1:m-1
  edges(cc,1:3)=[im*n,(im+1)*n,bright];
  cc++;
endfor

endfunction




