%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-11-26 21:02:20
 % @LastEditTime: 2022-01-03 23:19:44

 % get the bnid node idx for each cvt
 % NOTE: the bdofs must CONFORMS to phi
%--------------------------
function [bdofs_cell, dofid_cell, nbnodes_cell] = border_check(c_faces, c_nodes, c_vedges)

    cvt_num = size(c_faces, 1);
    bdofs_cell = cell(cvt_num, 1); % saves each bdofs
    dofid_cell = cell(cvt_num, 1); 

    nbnodes_cell = cell(cvt_num, 1); % saves the number of boundary dofs

%     figure;
% hold on;
    for ele = 1 : cvt_num        
        nodes = c_nodes{ele};
        vedges = c_vedges{ele};

        bnid = [];

        edge_num = size(vedges, 1);
        bnum = zeros(edge_num, 1);

        % NOTE: DO NOT using parfor! sequence related
        for i = 1 : edge_num
            n1 = vedges(i, 1:2);
            n2 = vedges(i, 3:4);

            bnid_each = [];
            for j = 1 : size(nodes, 1)
                p = nodes(j, :);
                
                flag = point_on_seg(n1, n2, p);
                if flag == 1
                    bnum(i) = bnum(i) + 1;
                    bnid_each = [bnid_each; j];
                end
            end
            bnid_new = rearrange_border_nodes(n1, nodes, bnid_each);     
            bnid = [bnid; bnid_new(1:end-1)]; % remove the adjoint corner node
%             bnid = [bnid; bnid_each(1:end-1)]; % remove the adjoint corner node
        end
        nbnodes_cell{ele} = bnum;

%         figure(1); clf;
%         scatter(nodes(bnid, 1), nodes(bnid, 2), 'filled'); hold on;
%         faces = c_faces{ele};
%         trimesh(faces, nodes(:,1), nodes(:,2));

        bdofs = [2*bnid-1, 2*bnid]';
        bdofs_cell{ele} = bdofs(:);

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

% rearrange as the distance to node n1
function bnid_new = rearrange_border_nodes(n1, nodes, bnid)

    nodes0 = nodes(bnid, :);

    dist = (nodes0(:,1) - n1(1)).^2 + (nodes0(:,2) - n1(2)).^2;
    nodes0(:,3) = dist;
    nodes0(:,4) = bnid;
    
    nodes1 = sortrows(nodes0, 3);
    bnid_new = nodes1(:, 4);
end

% segment = [n1, n2]
% if point p on the segment
function flag = point_on_seg(n1, n2, p)
    delta = 1e-2;
    if abs((p(1) - n1(1)) * (n2(2) - n1(2)) - (n2(1) - n1(1)) * (p(2) - n1(2))) < delta... % cross product
        && (min(n1(1), n2(1)) < p(1) || abs(min(n1(1), n2(1)) - p(1)) < delta)...
        && (p(1) < max(n1(1), n2(1)) || abs(p(1) - max(n1(1), n2(1))) < delta) ...
        && (min(n1(2), n2(2)) < p(2) || abs(min(n1(2), n2(2)) - p(2)) < delta) ...
        && (p(2) < max(n1(2), n2(2)) || abs(p(2) - max(n1(2), n2(2))) < delta)
        flag = 1;
    else
        flag = 0;
    end
end