%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-02-19 17:04:26
 % @LastEditTime: 2021-02-28 15:56:42
%--------------------------
function [iK, jK, edofMat] = prepare_assemble(nelx, nely, nelz)
    Num_node = (1+nely)*(1+nelx)*(1+nelz);
    nele = nelx*nely*nelz;
    
    nodenrs = reshape(1:Num_node,1+nely, 1+nelx, 1+nelz);
    edofVec = reshape(3*nodenrs(1:end-1,1:end-1,1:end-1)+1, nele, 1);
    edofMat = repmat(edofVec,1,24) + repmat([0 1 2 3*nely+[3 4 5 0 1 2] -3 -2 -1 ...
        3*(nelx+1)*(nely+1)+[0 1 2 3*nely+[3 4 5 0 1 2] -3 -2 -1]], nele, 1);

    iK = reshape(kron(edofMat,ones(24,1))',24*24*nele,1);
    jK = reshape(kron(edofMat,ones(1,24))',24*24*nele,1);
end