%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2022-01-24 18:50:51
 % @LastEditTime: 2022-04-28 11:22:51
%--------------------------
function [c_edofMat, c_nodes, cbns, c_nid, c_dupB] = assemble_cvt_free(c_polygon,)

    ncell = size(c_polygon, 1);
    c_nodes = cell(ncell, 1); % macro nodes of each cell
    c_nid   = cell(ncell, 1); % macro nodes indices of each cell and face
    c_dupB  = cell(ncell, 1); % macro nodes indices of each cell and face
    c_edofMat = cell(ncell, 1); % macro dofs of each cell

    nodes = [];
%     figure;

    % NOTE: no parfor, sequence dependent!
    for ele = 1 : ncell
        c_poly = c_polygon{ele};        
        [c_nid_f, c_dupB_f, nodesCell] = edof_one_cell(c_poly);

        c_nodes{ele} = nodesCell;    
        c_nid{ele}   = c_nid_f;    
        c_dupB{ele}  = c_dupB_f;

        nodes = [nodes; nodesCell];
    end
    cbns = uniquetol(nodes, 1e-6, "Byrows", true);
%     cbns = unique(nodes, "rows");
%     scatter3(nodes(:,1),nodes(:,2),nodes(:,3)); hold on;

    for ele = 1 : ncell

        nodesCell = c_nodes{ele};
        [~, idx] = ismembertol(nodesCell, cbns, 1e-5, 'Byrows', true); 
%         [~, idx] = ismember(nodesCell, cbns, 'rows'); 

        if sum(idx < 1) > 0
            error("here");
        end
        
        % NOTE: the sequence must conforms to Ke
        edof0 = [3*idx-2, 3*idx-1, 3*idx]';
        c_edofMat{ele} = edof0(:);
    end
end

function [c_nid, c_dupB, nodesCell] = edof_one_cell(c_poly)

    nface = size(c_poly, 1);
    c_nid = cell(nface, 1);
    c_dupB = cell(nface, 1);
    nodesCell = [];    

    % NOTE: no parfor, sequence dependent!
%     figure;
    for i = 1 : nface
        segsPerFace = c_poly{i}; % [n,6]

        [~, nodesFace] = edof_one_face(segsPerFace);

        % detect the duplicated nodes and remove them
        [c_dupB_node, nodesUni] = unique_control_points(nodesFace);

        nControlPnts = size(nodesUni, 1);

        % nid saves the idx of saved nodes in nodesCell
        if i==1
            nodesCell = nodesUni;
            nid = 1:nControlPnts;
        else
            % remove the duplicated nodes and add to nodesCell
            [uniIDX_face, nid] = deal_duplicated(nodesCell, nodesUni);
            nodesCell = [nodesCell; nodesUni(uniIDX_face,:)];  
        end             
        c_nid{i} = nid;
        c_dupB{i} = c_dupB_node;

%         scatter3( nodesUni(nid(:,2),1),  nodesUni(nid(:,2),2), nodesUni(nid(:,2),3), 'filled'); hold on;
%         x = segsPerFace(:,[1,4]);
%         y = segsPerFace(:,[2,5]);
%         z = segsPerFace(:,[3,6]);
%         plot3(x',y',z','LineWidth',2); hold on;  
    end
%     scatter3(nodesCell(:, 1), nodesCell(:, 2), nodesCell(:, 3), 'filled'); hold on;  
end

% TODO: HOW TO build cbns on surface
function [nedof, ref_nodes] = edof_one_face(segsPerFace)
    % detect end vertex polygon
    vertex = segsPerFace(:, 1:3);

    nseg = size(segsPerFace, 1);
    if nseg>7 && nseg~=10
        error("Side of Polygon > 7 or != 10 !");
% nseg
    end

    % s-patch
    switch nseg
    case 3
        % 200,101,002,011,020,110
        nvar = 6;
        ref_nodes = zeros(nvar, 3);
        ref_nodes([1,4,6],:) = vertex;
        ref_nodes(2,:) = mean(vertex([1,2],:));
        ref_nodes(3,:) = mean(vertex([1,3],:));
        ref_nodes(5,:) = mean(vertex([2,3],:));

    case 4
        % 2000,1100,1010,1001,0200,
        % 0110,0101,0020,0011,0002
        nvar = 10;
        ref_nodes = zeros(nvar, 3);
        ref_nodes([1,5,8,10],:) = vertex;
        ref_nodes(2,:) = mean(vertex([1,2],:));
        ref_nodes(3,:) = mean(vertex([1,3],:)); % 
        ref_nodes(4,:) = mean(vertex([1,4],:));
        ref_nodes(6,:) = mean(vertex([2,3],:));
        ref_nodes(7,:) = mean(vertex([2,4],:)); % 
        ref_nodes(9,:) = mean(vertex([3,4],:));

    case 5
        nvar = 15;
        ref_nodes = zeros(nvar, 3);
        ref_nodes(1,:) = vertex(1,:);
        ref_nodes(2,:) = mean(vertex([1,2],:));
        ref_nodes(3,:) = mean(vertex([1,3],:));
        ref_nodes(4,:) = mean(vertex([1,4],:));
        ref_nodes(5,:) = mean(vertex([1,5],:));
        ref_nodes(6,:) = vertex(2,:);
        ref_nodes(7,:) = mean(vertex([2,3],:));
        ref_nodes(8,:) = mean(vertex([2,4],:));
        ref_nodes(9,:) = mean(vertex([2,5],:));
        ref_nodes(10,:) = vertex(3,:);
        ref_nodes(11,:) = mean(vertex([3,4],:));
        ref_nodes(12,:) = mean(vertex([3,5],:));
        ref_nodes(13,:) = vertex(4,:);
        ref_nodes(14,:) = mean(vertex([4,5],:));
        ref_nodes(15,:) = vertex(5,:);

    case 6
        nvar = 21;
        ref_nodes = zeros(nvar, 3);
        ref_nodes(1,:) = vertex(1,:);
        ref_nodes(2,:) = mean(vertex([1,2],:));
        ref_nodes(3,:) = mean(vertex([1,3],:));
        ref_nodes(4,:) = mean(vertex([1,4],:)); % 
        ref_nodes(5,:) = mean(vertex([1,5],:));
        ref_nodes(6,:) = mean(vertex([1,6],:));
        ref_nodes(7,:) = vertex(2,:);
        ref_nodes(8,:) = mean(vertex([2,3],:));
        ref_nodes(9,:) = mean(vertex([2,4],:));
        ref_nodes(10,:) = mean(vertex([2,5],:)); % 
        ref_nodes(11,:) = mean(vertex([2,6],:));
        ref_nodes(12,:) = vertex(3,:);
        ref_nodes(13,:) = mean(vertex([3,4],:));
        ref_nodes(14,:) = mean(vertex([3,5],:));
        ref_nodes(15,:) = mean(vertex([3,6],:)); % 
        ref_nodes(16,:) = vertex(4,:);
        ref_nodes(17,:) = mean(vertex([4,5],:));
        ref_nodes(18,:) = mean(vertex([4,6],:));
        ref_nodes(19,:) = vertex(5,:);
        ref_nodes(20,:) = mean(vertex([5,6],:));
        ref_nodes(21,:) = vertex(6,:);

    case 7
        nvar = 28;        
        ref_nodes = zeros(nvar, 3);
        ref_nodes(1,:) = vertex(1,:);
        ref_nodes(2,:) = mean(vertex([1,2],:));
        ref_nodes(3,:) = mean(vertex([1,3],:));
        ref_nodes(4,:) = mean(vertex([1,4],:));
        ref_nodes(5,:) = mean(vertex([1,5],:));
        ref_nodes(6,:) = mean(vertex([1,6],:));
        ref_nodes(7,:) = mean(vertex([1,7],:));
        ref_nodes(8,:) = vertex(2,:);
        ref_nodes(9,:) = mean(vertex([2,3],:));
        ref_nodes(10,:) = mean(vertex([2,4],:));
        ref_nodes(11,:) = mean(vertex([2,5],:));
        ref_nodes(12,:) = mean(vertex([2,6],:));
        ref_nodes(13,:) = mean(vertex([2,7],:));
        ref_nodes(14,:) = vertex(3,:);
        ref_nodes(15,:) = mean(vertex([3,4],:));
        ref_nodes(16,:) = mean(vertex([3,5],:));
        ref_nodes(17,:) = mean(vertex([3,6],:));
        ref_nodes(18,:) = mean(vertex([3,7],:));
        ref_nodes(19,:) = vertex(4,:);
        ref_nodes(20,:) = mean(vertex([4,5],:));
        ref_nodes(21,:) = mean(vertex([4,6],:));
        ref_nodes(22,:) = mean(vertex([4,7],:));
        ref_nodes(23,:) = vertex(5,:);
        ref_nodes(24,:) = mean(vertex([5,6],:));
        ref_nodes(25,:) = mean(vertex([5,7],:));
        ref_nodes(26,:) = vertex(6,:);
        ref_nodes(27,:) = mean(vertex([6,7],:));
        ref_nodes(28,:) = vertex(7,:);
        
    case 10
        nvar = 55;        
        ref_nodes = zeros(55, 3);
        idx = 1;
        for i = 1:10
            for j = i:10
                ref_nodes(idx, :) = mean(vertex([i,j],:));
                idx = idx + 1;
            end
        end
    end
    nedof = nvar * 3;
end

function [uniIDX_local, nid] = deal_duplicated(nodesGlobal, nodesLocal)

    tol = 1e-5;

    nid = []; 
    uniIDX_local = [];
    k = size(nodesGlobal, 1) + 1;

    for i = 1 : size(nodesLocal, 1)
        [idx, ~] = find(abs(nodesGlobal(:, 1)-nodesLocal(i, 1)) < tol & ...
                        abs(nodesGlobal(:, 2)-nodesLocal(i, 2)) < tol & ...
                        abs(nodesGlobal(:, 3)-nodesLocal(i, 3)) < tol);

        if isempty(idx)
            uniIDX_local = [uniIDX_local; i];
            nid = [nid; k];
            k = k+1;
        else
            if length(idx) > 1
                error("here");
            end
            nid = [nid; idx];
        end
    end
end

function [dupID, nodesUni] = unique_control_points(nodesFace)

    % remove the duplicated nodes
    nodesUni0 = uniquetol(nodesFace, 1e-5, "Byrows", true);

    % reorder the unique nodes as original order
    [nodesUni, uniID, ~] = intersect(nodesFace, nodesUni0, 'rows', 'stable');

    % find the duplicated idx in nodesFace
    tol   = 1e-5;
    dupID = cell(size(nodesUni, 1), 1);

    for i = 1 : size(nodesUni, 1)
        [idx, ~] = find(abs(nodesFace(:, 1)-nodesUni(i, 1)) < tol & ...
                        abs(nodesFace(:, 2)-nodesUni(i, 2)) < tol & ...
                        abs(nodesFace(:, 3)-nodesUni(i, 3)) < tol);
        dupID{i} = idx; % length(idx) may be 3
    end
end