%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-01-25 14:11:35
 % @LastEditTime: 2021-01-28 21:27:38

 % given microx, parameterize B of 4 boundarys based on coordinates
 % Both the order of interpolated micro boundary and macro variables are:
 % [bottom -> right -> up -> left]
%--------------------------

function [B_boundary, seg_num] = parametric_B(microx, addNum)

    bdofs_num = 2*(microx+1); 
    border  = round(linspace(1, microx+1, addNum+2)); 
    seg_num = length(border)-1; % segment num divided by added macro dofs
    
    % macro dofs num of ONE boundary, e.g. u_c1,p1,p2,u_b,p3,p4,uc2
    % the dofs of u_b belong to seg of [u_c1,p1,p2,u_b] and [u_b,p3,p4,uc2]
    % seg_num*6 + 2 could remove duplicated dofs 
    % if addNum = 0, seg_num = 1, var_num_1b = 8
    var_num_1b = seg_num*6 + 2;   
    var_num_4b = 4 * (var_num_1b - 2);
    B_1b = zeros(bdofs_num, var_num_1b);
    
    % for every segment of one boundary, B maybe different
    ir = 1; ic = 1;
    for i = 1 : seg_num-1
        nodes_num = border(i+1) - border(i) + 1;
        B = parametric_B_one_boundary(nodes_num); % [2,8,x], x is nodes num of every seg

        % reshape B from [2, 8, x+1] to [2(x+1), 8]
        B1 = permute(B, [2,1,3]); % [8, 2, x+1]
        B2 = reshape(B1, 8, [])'; % [2(x+1), 8]

        % remove duplicated nodes caused by u_b
        % here 8 is caused by P_num = 2
        B_1b(ir:ir + nodes_num*2-2-1, ic:ic+8-1) = B2(1:end-2, :);

        ir = ir + nodes_num*2-2;
        ic = ic + 8-2;
    end   
    % assemble B_1b with the last seg  
    i = seg_num;
    nodes_num = border(i+1) - border(i) + 1;
    B = parametric_B_one_boundary(nodes_num); % [2,8,x], x is nodes num of every seg
    B1 = permute(B, [2,1,3]); % [8, 2, x+1]
    B2 = reshape(B1, 8, [])'; % [2(x+1), 8]

    % The last seg doesn't have duplicated nodes
    B_1b(ir:ir + nodes_num*2-1, ic:ic+8-1) = B2;

    % 4 boundarys, NOTE: the nodes sequence  
    B_bottom = [B_1b, zeros(bdofs_num, var_num_4b-var_num_1b)];
    B_right  = [zeros(bdofs_num, var_num_1b-2), B_1b, zeros(bdofs_num, var_num_1b*2-6)];
    B_up     = [zeros(bdofs_num, var_num_1b*2-4), B_1b, zeros(bdofs_num, var_num_1b-4)];
    B_left   = [B_1b(:, end-1:end), zeros(bdofs_num, var_num_4b-var_num_1b), B_1b(:, 1:end-2)];    
    
    % NOTE: remove duplicated nodes
    B_boundary = [B_bottom; B_right(3:end, :); B_up(3:end, :); B_left(3:end-2, :)];
end

% nodes_num is the num of nodes on this boundary
function B = parametric_B_one_boundary(nodes_num)
    
    B = zeros(2, 8, nodes_num);

    for i = 0 : nodes_num-1
        t = i / (nodes_num-1); % project to [0,1]
        B(:, :, i+1) = Bezier_matrix(t);
    end
end

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

% test B
% addNum = 2;
% [~, ~, edofMat_ma1, ~] = macro_edofMat_bezier(nelx, nely, addNum);
% [ref_nodes, ~, ~] = multi_generate_mesh(lx, nelx, nely, addNum);
% U = reshape(ref_nodes', [], 1);

% addNum = microx - 1;
% [~, ~, edofMat_ma2, ~] = macro_edofMat_bezier(nelx, nely, microx-1);
% alldofs = 2*(nelx+1)*(nely+1) + 2 * addNum*nely*(nelx+1) + 2 * addNum*nelx*(nely+1);
% U2 = zeros(alldofs, 1);
% for ele = 1:eleNum_ma
%     edof = edofMat_ma1(ele, :);
%     ub = bezier_B * U(edof);

%     edof = edofMat_ma2(ele, :);
%     U2(edof) = ub;
% end
% deform_nodes = reshape(U2, 2, [])';
% figure; clf; scatter(deform_nodes(:,1), deform_nodes(:,2), 'b'); axis equal;
% hold on; scatter(ref_nodes(:,1), ref_nodes(:,2), 'k', 'filled'); axis equal;