%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [q0,edof] = mesh_sqr(xbox,ybox,nx,ny)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b = xbox(2)-xbox(1);
h = ybox(2)-ybox(1);

nd = 1;

for xx = 0:1:nx
  for yy = 0:1:ny
    x0(nd,1) = xx*b/nx;
    y0(nd,1) = yy*h/ny;
    nd = nd +1;
  end
end  

[size0,sizen]=size(x0);

q0(1:2:2*size0,1) = x0(:);
q0(2:2:2*size0,1) = y0(:);

[ndof,sizen]=size(q0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% regular quads %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ie = 1;
for xx = 1:1:nx
  for yy = 1:1:ny
    nn = (ny+1)*(xx-1)+yy;
    el0(ie,1) = nn;
    el0(ie,2) = nn+ny+1;
    el0(ie,3) = nn+ny+2;
    el0(ie,4) = nn+1;
    ie = ie +1;
  end
end  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% edof ... global incidence matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nel,sizen]=size(el0);

for i=1:nel
 edof(i,1)=i;
 for j=1:4
  ic1=2*j; ic2=2*j+1; 
  edof(i,ic1)=el0(i,j)*2-1;
  edof(i,ic2)=el0(i,j)*2;
 end
end
size(edof)
edof
disp(['*** generated ',num2str(nel),' *** elems with ',num2str(ndof),' dofs ***']);      

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
