function [node, element]=  c2d9(a, ndiv)
% author: Ji Wan, Wuhan University, Email: wanji@whu.edu.cn
% date: 2021\06\25
npx = 2*ndiv(1)+1;
npy = 2*ndiv(2)+1;
[x, y] = meshgrid(linspace(0, a(1), npx), ...
    linspace(0, a(2), npy));
%% eight nodes
[p1, q1] = meshgrid(1:2:npx-1,1:2:npy-1);
pos1 = sub2ind([npx, npy], p1, q1);
[p1, q1] = meshgrid(3:2:npx,1:2:npy-1);
pos2 = sub2ind([npx, npy], p1, q1);
[p1, q1] = meshgrid(3:2:npx,3:2:npy);
pos3 = sub2ind([npx, npy], p1, q1);
[p1, q1] = meshgrid(1:2:npx-1,3:2:npy);
pos4 = sub2ind([npx, npy], p1, q1);
[p1, q1] = meshgrid(2:2:npx-1,1:2:npy-1);
pos5 = sub2ind([npx, npy], p1, q1);
[p1, q1] = meshgrid(3:2:npx,2:2:npy);
pos6 = sub2ind([npx, npy], p1, q1);
[p1, q1] = meshgrid(2:2:npx, 3:2:npy);
pos7 = sub2ind([npx, npy], p1, q1);
[p1, q1] = meshgrid(1:2:npx-1, 2:2:npy);
pos8 = sub2ind([npx, npy], p1, q1);
[p1, q1] = meshgrid(2:2:npx, 2:2:npy);
pos9 = sub2ind([npx, npy], p1, q1);
element = [reshape(pos1, numel(pos1), 1),... 
    reshape(pos2, numel(pos2), 1),...
    reshape(pos3, numel(pos3), 1),...
    reshape(pos4, numel(pos4), 1),...
    reshape(pos5, numel(pos5), 1),...
    reshape(pos6, numel(pos6), 1),...
    reshape(pos7, numel(pos7), 1),...
    reshape(pos8, numel(pos8), 1),...
    reshape(pos9, numel(pos9), 1)];
node = [reshape(x',numel(x),1), reshape(y', numel(y),1)];
end