%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-03-01 19:44:43
 % @LastEditTime: 2021-03-05 21:14:43
%--------------------------
function N = trilinear_interpolation_border(microx, bdofs)
    Num_node = (1+microx)^3;
    nodenrs = reshape(1:Num_node, 1+microx, 1+microx, 1+microx);

    bdofs_num = length(bdofs);
    N = zeros(bdofs_num, 24);

    bnodes = bdofs(3:3:end)/3;
    for i = 1 : bdofs_num/3
        % compute coor of this border nodes
        [row, col] = find(nodenrs == bnodes(i));
        [ix, iz] = row_col_idx(microx+1, col);

        cy = row-1; % [0, microx]
        cx = ix-1;
        cz = iz-1;
        Nb = trilinear_ele(microx, cx, microx-cy, cz);

        dofid = [3*i-2, 3*i-1, 3*i];
        N(dofid, :) = Nb;
    end
end

 %   y              4-----3
 %   |             /|    /|
 %   0---x        8-----7 2   
 %  /             |     |/
 % z              5-----6
function Nb = trilinear_ele(lx, x, y, z)
    microx = 1;
    Num_node = (1+microx)^3;
    nodenrs = reshape(1:Num_node, 1+microx, 1+microx, 1+microx);    
    
    N(1) = (1-x/lx) * (1-y/lx) * (1-z/lx);
    N(2) = x/lx * (1-y/lx) * (1-z/lx);
    N(3) = x/lx * y/lx * (1-z/lx);
    N(4) = (1-x/lx) * y/lx * (1-z/lx);
    N(5) = (1-x/lx) * (1-y/lx) * z/lx;
    N(6) = x/lx * (1-y/lx) * z/lx;
    N(7) = x/lx * y/lx * z/lx;
    N(8) = (1-x/lx) * y/lx * z/lx;

    % reorder N to macro-dofs' order
    nodes = N;
    nodes(2, :) = [nodenrs(2,1,1); nodenrs(2,2,1); nodenrs(1,2,1); nodenrs(1,1,1); 
                   nodenrs(2,1,2); nodenrs(2,2,2); nodenrs(1,2,2); nodenrs(1,1,2);]';
    nodes = sortrows(nodes', 2);
    N = nodes(:, 1);

    Nb = [N(1), 0, 0, N(2), 0, 0, N(3), 0, 0, N(4), 0, 0, N(5), 0, 0, N(6), 0, 0, N(7), 0, 0, N(8), 0, 0;
            0, N(1), 0, 0, N(2), 0, 0, N(3), 0, 0, N(4), 0, 0, N(5), 0, 0, N(6), 0, 0, N(7), 0, 0, N(8), 0;
            0, 0, N(1), 0, 0, N(2), 0, 0, N(3), 0, 0, N(4), 0, 0, N(5), 0, 0, N(6), 0, 0, N(7), 0, 0, N(8)];
end

% % test
% u0 = reshape(ref_nodes', [], 1);

% num_b = microx-1;
% ref_nodes = mesh3D(nelx, nely, nelz, microx, num_b, lx, opt_scale);
% u1 = reshape(ref_nodes', [], 1);
% u2 = zeros(size(u1));
% [iK_ma, jK_ma, edofMat2] = assemble_bezier(nelx, nely, nelz, num_b);

% dofs = sort([corners; add_dofs]);
% [corner_id,~] = find(dofs==corners');

% figure;
% for ele = 1:nele_ma
%     edof = edofMat_ma(ele, :);
%     ue = u0(edof);
%     uc = ue(corner_id);

%     ub = reshape(N * uc, 3, [])';
%     scatter3(ub(:,1), ub(:,2), ub(:,3),'filled'); axis equal; hold on;
%     uc = reshape(uc, 3, [])';
%     scatter3(uc(:,1), uc(:,2), uc(:,3),'filled'); axis equal; hold on;
%     ue = reshape(ue, 3, [])';
%     scatter3(ue(:,1), ue(:,2), ue(:,3),'filled'); axis equal; hold on;
% end