function def_nodes = extract_def_nodes_macro(U, nelx, nely, nelz, lx, num_b, edofMat_control, opt_load_edofmat)

    [corners, boundarys, ~] = micro_boundary2(num_b+1);
    surfaces = [corners(:); boundarys(:)];
    dofs = sort(surfaces);
    [dofid,~,~] = find(dofs==corners');
    
    ref_nodes = mesh3D(nelx, nely, nelz, 1, 0, lx, 1);

    if nargin > 7
        if opt_load_edofmat==1
            load dizhi_edofMat0
        else
            [~, ~, edofMat_macro] = assemble_bezier(nelx, nely, nelz, 0);
            save dizhi_edofMat0 edofMat_macro
        end
    else
        [~, ~, edofMat_macro] = assemble_bezier(nelx, nely, nelz, 0);
    end

    udef = zeros(3*(nelx+1)*(nely+1)*(nelz+1), 1);
    for ele = 1 : nelx*nely*nelz
        edof = edofMat_control(ele, :);
        ue = U(edof);

        edof = edofMat_macro(ele, :);
        udef(edof) = ue(dofid); 
    end
    def_nodes = ref_nodes + reshape(udef, 3, [])';
%     figure; clf; scatter3(def_nodes(:,1), def_nodes(:,2), def_nodes(:,3), 'filled'); axis equal; 
end

function [cdofs, bdofs, idofs] = micro_boundary2(microx)
    Num_node = (1+microx)^3;
    nodenrs = reshape(1:Num_node, 1+microx, 1+microx, 1+microx);

    corners = nodenrs([1,end],[1,end],[1,end]);
        
    % extract surface nodes
    inner = nodenrs(2:end-1, 2:end-1, 2:end-1);
    bb = setdiff(nodenrs(:), inner(:));

    boundarys = setdiff(bb(:), corners(:));
    
    bdofs = nodes2dofs(boundarys);
    idofs = nodes2dofs(inner);
    cdofs = nodes2dofs(corners);
end