function [nodes,edges,elemNodes,elemEdges,leftNodes,topNodes,rightNodes,bottomNodes,leftEdges,topEdges,rightEdges,bottomEdges] = rectangulateRectangle(a,b,hx,hz)

% [nodes,edges,elemNodes,elemEdges,leftNodes,topNodes,rightNodes,bottomNodes,leftEdges,topEdges,rightEdges,bottomEdges] = rectangulateRectangle(a,b,hx,hz)

m1 = ceil(a/hx);
n1 = m1+1;
x1 = linspace(0,a,n1);

m2 = ceil(b/hz);
n2 = m2+1;
x2 = linspace(0,b,n2);

[X1,X2] = meshgrid(x1,x2);
nodes = [X1(:)';X2(:)'];
n = n1*n2;

me = (m1+1)*m2+(m2+1)*m1;
edges = zeros(2,me);

m = m1*m2;
elemNodes = zeros(4,m);
elemEdges = zeros(4,m);

idx = 1;
idxEdges = 1;
for i=1:m1
    for j=1:m2
        edges(:,idxEdges+[0,1]) = [ [(i-1)*n2+j;i*n2+j] , [(i-1)*n2+j;(i-1)*n2+j+1] ];
        elemNodes(:,idx) = [(i-1)*n2+j;i*n2+j;i*n2+j+1;(i-1)*n2+j+1];
        if i<m1 && j<=m2
            elemEdges(:,idx) = [idxEdges;idxEdges+2;idxEdges+1;idxEdges+2*m2+2];
        elseif i==m1
            elemEdges(:,idx) = [idxEdges;idxEdges+2;idxEdges+1;me-m2+j];
        end
        idxEdges = idxEdges+2;
        idx = idx+1;
    end
    edges(:,idxEdges) = [i*n2;(i+1)*n2];
    idxEdges = idxEdges+1;
end
for i=1:m2
    edges(:,idxEdges) = [n-n2+i;n-n2+i+1];
    idxEdges = idxEdges+1;
end

leftNodes = 1:n2;
topNodes = n2:n2:n;
rightNodes = n-n2+1:n;
bottomNodes = 1:n2:n-n2+1;

leftEdges = -[2*m2:-2:2];
rightEdges = me-m2+1:me;
bottomEdges = 1:2*m2+1:me-m2;
topEdges = -[me-m2:-(2*m2+1):2*m2+1];
