%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-01-25 16:54:13
 % @LastEditTime: 2021-01-25 16:57:26
%--------------------------
% order: col by col
function [iK, jK, multi_edofMat, EleNodesID] = macro_edofMat_bezier(nelx, nely, addNum)

    dofs_num = 4*2 + addNum*4*2;
    multi_edofMat = zeros(nelx*nely, dofs_num);

    for i = 1:nelx
        for j = 1:nely
            ele = j + (i-1)*nely;
            multi_edof = interpolated_edof(i, j, addNum, nely);
            multi_edofMat(ele, :) = multi_edof;            
        end
    end

    iK = reshape(kron(multi_edofMat,ones(dofs_num,1)).', dofs_num^2*nelx*nely, 1);
    jK = reshape(kron(multi_edofMat,ones(1,dofs_num)).', dofs_num^2*nelx*nely, 1);

    EleNodesID = multi_edofMat(:, 2:2:end)./2; 
end

% every macro ele: down -> right -> up -> left
% microx: controllable nodes num, must > 0
% this ele idx: [ely, elx]
%      n-n-n-n-n-n-n
% e.g. n     n     n
%      n     n     n
%      n-n-n-n-n-n-n
% ynum_v h+h v h+h v
function edof = interpolated_edof(elx, ely, addNum, nely)
    
    % 1. get nodes idx
    
    ynum_v = nely+1 + addNum*nely; % node num in vertical edge in y-dir
    ynum_h = (nely+1) * addNum; % node num across horizontal edge in y-dir among 1 macro-ele range    

    % idx: last + current
    left_up_node  = (ynum_v + ynum_h)*(elx-1) + (ely-1)*(1 + addNum) + 1;
    right_up_node = (ynum_v + ynum_h)* elx    + (ely-1)*(1 + addNum) + 1;

    corners = [left_up_node+addNum+1, right_up_node+addNum+1, right_up_node, left_up_node];    

    first_up_node = ynum_v*elx + ynum_h*(elx-1) + ely;
    if addNum<2
        last_up_node = first_up_node;
    else
        last_up_node = ynum_v*elx + ynum_h*(elx-1) + (nely+1) * (addNum-1) + ely;
    end
    
    down = first_up_node+1 : nely+1 : last_up_node+1;
    up   = first_up_node   : nely+1 : last_up_node;    

    right = right_up_node+1 : right_up_node+addNum;
    left  = left_up_node +1 : left_up_node +addNum;
    
    % 2. post-pocessing & nodes idx -> dofs
    up    = flip(up);
    right = flip(right);
    
    up_d      = node2dof(up);
    down_d    = node2dof(down);
    left_d    = node2dof(left);
    right_d   = node2dof(right);
    corners_d = node2dof(corners);

    edof = [corners_d([1,2]), down_d, corners_d([3,4]), right_d, corners_d([5,6]), up_d, corners_d([7,8]), left_d];
end

% node is 
function dof = node2dof(node)
    idx = [2*node(:)-1, 2*node(:)]';
    dof = idx(:)';
end