%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-11-26 21:02:20
 % @LastEditTime: 2022-04-28 10:47:50

 % for each cvt and each surface: get the boundary node idx
 % NOTE: the bdofs must CONFORMS to phi
 % CHECK: bdofs order, do we need rearrange?
 % paras:
 % c_polygon: all segs for each face and each cell
%--------------------------
function [c_bnid, c_dofid, c_bnidSaved, c_bnidSavedID, nbnidFine] = border_check(c_node, c_polygon)

    ncell = size(c_node, 1);
    c_dofid = cell(ncell, 1);
    c_bnid  = cell(ncell, 1);
    c_bnidSavedID = cell(ncell, 1);
    c_bnidSaved   = cell(ncell, 1);
    nbnidFine = zeros(ncell, 1);

%     figure;
    for ele = 1 : ncell
        nodes = c_node{ele};
        polys = c_polygon{ele};    

        [c_bnid_face, c_bnidSaved_face, c_bnidSavedID_face, bnid] = boundary_nodes(nodes, polys);

        c_bnid{ele} = c_bnid_face;
        c_bnidSavedID{ele} = c_bnidSavedID_face;
        c_bnidSaved{ele}   = c_bnidSaved_face;

        nbnidFine(ele) = length(bnid);
%         bnid = sort(bnid); % important, phi-B is order as this
        bdofs = [3*bnid-2, 3*bnid-1, 3*bnid]';

        alldofs = 1 : size(nodes, 1) * 3;
        idofs   = setdiff(alldofs, bdofs(:));
        newdofs = [bdofs(:)', idofs];
        dofid   = reorder_dofs(newdofs); 
        c_dofid{ele} = dofid;        
    end
end

function [c_bnid, c_bnidSaved, c_bnidSavedID, bnid] = boundary_nodes(nodes, c_poly)

    tol = 1e-3;

    nface = size(c_poly, 1);
    npnt  = size(nodes, 1);
                             
    % bnid for each face 
    c_bnid = cell(nface, 1); 
    % bnid for each face, removed duplicated node with previous
    c_bnidSaved = cell(nface, 1); 
    % index of bnidSaved, for phi-B
    c_bnidSavedID = cell(nface, 1);
    bnid = [];

%     figure; scatter3(nodes(:,1), nodes(:,2), nodes(:,3),'b'); hold on;
    % NOTE: DO NOT using parfor! sequence related
    for i = 1 : nface
        segs = c_poly{i};

        p1 = segs(1, 1:3);
        v1 = segs(1, [4:6]) - segs(1, [1:3]);
        v2 = segs(2, [4:6]) - segs(2, [1:3]);

        % TODO: change to surface

        plane = repmat([p1, v1, v2], npnt, 1);
        projNode = projPointOnPlane(nodes, plane);

        dist = projNode - nodes;
        [idx, ~] = find(abs(dist(:, 1)) < tol & ...
                        abs(dist(:, 2)) < tol & ...
                        abs(dist(:, 3)) < tol);

        c_bnid{i} = idx; % bnid of this face

        [bnidSaved, ia] = setdiff(idx, bnid, 'stable');
        bnid = [bnid; bnidSaved(:)];

        % id of bnidSaved of this face
%         [~,bnidSavedIDX,~] = intersect(idx, bnidSaved, 'rows', 'stable');

        c_bnidSavedID{i} = ia;
        c_bnidSaved{i}   = length(bnid) - length(bnidSaved) + 1 : length(bnid);

%         x = segs(:,[1,4]);
%         y = segs(:,[2,5]);
%         z = segs(:,[3,6]);
%         plot3(x',y',z','LineWidth',2); hold on;  
    end     
%     scatter3(nodes(bnid,1), nodes(bnid,2), nodes(bnid,3), 'filled', 'r'); hold on;
end