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

    % extract surface nodes
    inner = nodenrs(2:end-1, 2:end-1, 2:end-1);
    surface_nodes = setdiff(nodenrs(:), inner(:));
%     corners = nodenrs([1,end],[1,end],[1,end]);
%     surface_nodes = setdiff(surface_nodes(:), corners(:));
    dofs_num = length(surface_nodes)*3;
    N = zeros(dofs_num, 24);

    for i = 1 : microx+1
        for j = 1 : microx+1
            for k = [1, microx+1]                
                Nb = trilinear_ele(microx, i-1, j-1, k-1);
                nid = nodenrs(j,i,k);
                dof = get_dofs_id(surface_nodes, nid);
                N(dof, :) = Nb;
            end
        end
    end
    for i = [1, microx+1]
        for j = 1 : microx+1
            for k = 1 : microx+1
                Nb = trilinear_ele(microx, i-1, j-1, k-1);
                nid = nodenrs(j,i,k);
                dof = get_dofs_id(surface_nodes, nid);
                N(dof, :) = Nb;
            end
        end
    end
    for i = 1 : microx+1
        for j = [1, microx+1]
            for k = 1 : microx+1
                Nb = trilinear_ele(microx, i-1, j-1, k-1);
                nid = nodenrs(j,i,k);
                dof = get_dofs_id(surface_nodes, nid);
                N(dof, :) = Nb;
            end
        end
    end
end

function dofs_id = get_dofs_id(all_nodes, nid)
    [nid1, ~] = find(all_nodes==nid);
    dofs_id = [3*nid1(:)-2, 3*nid1(:)-1, 3*nid1(:)]';
    dofs_id = dofs_id(:);
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(1,1,1);
    %                nodenrs(2,1,2); nodenrs(2,2,2); nodenrs(1,2,2); nodenrs(1,1,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
% [boundarys, inner] = micro_boundary(microx);
% 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);

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

%     ub = [ N] * ue;
%     ub(:,2) = boundarys;
%     ub = sortrows(ub,2);
%     edof = edofMat2(ele, :);
%     u2(edof) = ub(:,1);
% end
% nodes = reshape(u2, 3, [])';
% figure; scatter3(nodes(:,1), nodes(:,2), nodes(:,3)); axis equal;