%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-11-27 15:21:49
 % @LastEditTime: 2021-12-06 20:12:34

 % could be arbitrary version without considering the length?
 % NOTE: the bdofs must CONFORMS to phi
 % bdofs order of phi: boundary dofs on each edge
 % Three situations:
 % 1. fine boundary nodes >= 2 on this edge, phi as usual
 % 2. fine boundary node   = 1 on this edge, phi with 1 contro node
 % 3. fine boundary node   = 0 on this edge, phi = eye(4)
%--------------------------
function phi_cell = parametric_B_cvt(nbnodes_cell, nbridge_each)

    ncvt = size(nbnodes_cell, 1);
    phi_cell = cell(ncvt,1);
%     figure;

    B2 = simple_quad_Bezier(3);
    
    for ele = 1 : ncvt
        nbele = nbnodes_cell{ele} - 1;
        nedges = size(nbele, 1);

        ncbn = 1;
        nbdofs = 1;

        standard_nedges = length(find(nbele>2));
        ncbn_1v = standard_nedges * ((nbridge_each+1) * 2 + nbridge_each) + nedges;
        ncbn_1v = ncbn_1v + sum(nbele(nbele <=2) - 1); % for situations 2 & 3
        phi_1v = zeros(sum(nbele)*2, ncbn_1v * 2);

        for edge = 1 : nedges
            nb_edge = nbele(edge);

            % three situations
            if nb_edge > 2
                % 1. fine boundary nodes >= 2 on this edge, phi as usual
                [phi_1b, seg_num] = bezier_one_edge(nb_edge, nbridge_each);
                ncbn_1b = seg_num * 2 + nbridge_each + 2; 
            elseif nb_edge == 2
                % 2. fine boundary node   = 1 on this edge, phi with 1 contro node
                % linear interpolation
                phi_1b = B2;
                ncbn_1b = 3; 
            else
                % 3. fine boundary node   = 0 on this edge, phi = ones(2,2)
                phi_1b = eye(4);
                ncbn_1b = 2; 
            end
            
            % to remove the adjoint corner_node
            if edge == nedges
                phi_1v(nbdofs : nbdofs + nb_edge*2 -1, ncbn : ncbn + (ncbn_1b-1)*2 - 1) = phi_1b(1:end-2, 1:end-2);
                phi_1v(nbdofs : nbdofs + nb_edge*2 -1, 1:2) = phi_1b(1:end-2, end-1:end);
            else
                phi_1v(nbdofs : nbdofs + (nb_edge+1)*2 -1, ncbn : ncbn + ncbn_1b*2 - 1) = phi_1b;
            end

            % to remove the adjoint corner_node
            ncbn = ncbn + (ncbn_1b-1) * 2;
            
            nbdofs = nbdofs + nb_edge * 2;
        end
        phi_cell{ele} = phi_1v;
    end
end

function B2 = simple_quad_Bezier(nodes_num)
    B = zeros(2, 6, nodes_num);

    for i = 0 : nodes_num-1
        t = i / (nodes_num-1); % project to [0,1]
        B(:, :, i+1) = Bezier_matrix(t);
    end
        
    B1 = permute(B, [2,1,3]); % [6, 2, x+1]
    B2 = reshape(B1, 6, [])'; % [2(x+1), 6]
end

function B = Bezier_matrix(t)
    b = [1, t, t^2] * [1, 0,0;
                        -2,2,0;
                        1,-2,1]; % 1,3
    I = eye(2, 2);
    B = kron(b, I); % 2,6
end