function [node, element]=  c2d4(a, ndiv)
% author: Ji Wan, Wuhan University, Email: wanji@whu.edu.cn
% date: 2021-06-25
[x, y] = meshgrid(linspace(0, a(1), ndiv(1)+1), ...
    linspace(0, a(2), ndiv(2)+1));
%% four nodes
[p1, q1] = meshgrid(1:ndiv(1),1:ndiv(2));
pos1 = sub2ind([ndiv(1)+1, ndiv(2)+1], p1, q1);
[p2, q2] = meshgrid(2:ndiv(1)+1,1:ndiv(2));
pos2 = sub2ind([ndiv(1)+1, ndiv(2)+1], p2, q2);
[p3, q3] = meshgrid(2:ndiv(1)+1,2:ndiv(2)+1);
pos3 = sub2ind([ndiv(1)+1, ndiv(2)+1], p3, q3);
[p4, q4] = meshgrid(1:ndiv(1),2:ndiv(2)+1);
pos4 = sub2ind([ndiv(1)+1, ndiv(2)+1], p4, q4);
element = [reshape(pos1, numel(pos1), 1),... 
    reshape(pos2, numel(pos2), 1),...
    reshape(pos3, numel(pos3), 1),...
    reshape(pos4, numel(pos4), 1)];
node = [reshape(x',numel(x),1), reshape(y', numel(y),1)];
end