%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-10-29 15:04:18
 % @LastEditTime: 2021-02-18 22:22:38
%--------------------------
% dofs order: col by col
% outer order: 4_coners, bdofs:[bottom, right, up, left]
% idofs: border-dofs should be intepolated
% border = [outer, idofs, fdofs]
% addNum could be zero
function [coners, bdofs, idofs, inner, addNum] = multi_microBC(nelx, addNum, opt_itp)  

    nodenrs = reshape(1:(1 + nelx)^2, 1 + nelx, 1 + nelx);    
    alldofs = 1:2 * (nelx + 1)^2;
    
    if addNum<3
        border = round(linspace(1, nelx+1, addNum+2));
    else
        tmp = round(linspace(2, nelx, addNum));
        border = [1, tmp, nelx+1];
    end

    add_indices = unique(border(2:end-1)); % NOTE: remove duplicate elements
    addNum = length(add_indices); % get new addNum, and return

    % four corner points
    n1 = [nodenrs(end, [1, end]), nodenrs(1, [end, 1])];
    dcorner = reshape([2 * n1 - 1; 2 * n1], 1, 8); 

    % bottom boundaries
    nbottom = nodenrs(end, add_indices);
    dbottom = reshape([2 * nbottom - 1; 2 * nbottom], 1, 2*addNum); 

    nbottom_full = nodenrs(end, 2:end - 1);
    dbottom_full = reshape([2 * nbottom_full - 1; 2 * nbottom_full], 1, []); 
    dbottom_b = setdiff(dbottom_full, dbottom);

    % right boundaries, NOTE: REVERSE!
    nright = flip(nodenrs(add_indices, end))';
    dright = reshape([2 * nright - 1; 2 * nright], 1, 2*addNum);

    nright_full = nodenrs(2:end - 1, end);
    dright_full = reshape([2 * nright_full - 1; 2 * nright_full], 1, []); 
    dright_b = setdiff(dright_full, dright);    
    tmp = fliplr(reshape(dright_b,2,[])); % reverse
    dright_b = tmp(:)';
    
    % up, NOTE: REVERSE!
    nup = flip(nodenrs(1, add_indices));
    dup = reshape([2 * nup - 1; 2 * nup], 1, 2*addNum); 

    nup_full = nodenrs(1, 2:end - 1);
    dup_full = reshape([2 * nup_full - 1; 2 * nup_full], 1, []); 
    dup_b = setdiff(dup_full, dup);    
    tmp = fliplr(reshape(dup_b,2,[])); % reverse
    dup_b = tmp(:)';

    % left
    nleft = nodenrs(add_indices, 1)';
    dleft = reshape([2 * nleft - 1; 2 * nleft], 1, 2*addNum); 

    nleft_full = nodenrs(2:end - 1, 1);
    dleft_full = reshape([2 * nleft_full - 1; 2 * nleft_full], 1, []); 
    dleft_b = setdiff(dleft_full, dleft);

    coners = dcorner(:);
    bdofs = [dbottom, dright, dup, dleft]';
    outer = [coners; bdofs]; 
        
    if opt_itp
        idofs = [dbottom_b, dright_b, dup_b, dleft_b]'; 
        inner = setdiff(alldofs, [outer; idofs])';
    else
        idofs = [];
        inner = setdiff(alldofs, outer)';
    end
end