%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-03-31 20:20:38
 % @LastEditTime: 2021-08-16 14:57:10
%--------------------------
function [regulated_matrix, dFdx, B, NH, sK] = precomputation(eleNum_ma, microx, rho, bezier_B, ke0, ke1, iK_mi, jK_mi, ...
    bdofs_num, dofid, edofMat_mi, dx, nelx, nely, num_b)

    eleNum_mi = microx^2;
    regulated_matrix  = cell(eleNum_ma, 1);
    B    = cell(eleNum_ma, 1);
    dFdx = cell(4, eleNum_ma);
%     I    = speye(bdofs_num);
    
    bezier_B = sparse(bezier_B);
    order1 = reshape(edofMat_mi', [], 1);
    
    vdofs_num = size(bezier_B, 2);

    % NOTE: Nh is ordered as 
    % 2 4
    % 1 3
    [dNh_ele, ~, Nh]   = shape_gradient_ele(dx/2, dx/2, 4);     
    
    lx = dx * microx;
    NH = cell(eleNum_ma, 1); 
    
    sK = zeros(vdofs_num^2, eleNum_ma);
        
    for ele = 1 : eleNum_ma          
%         tic
        material = rho{ele};
%         material(:) = 1;

        % Assemble the micro global stiffness as the order of [dofid]
        % i.e. [4-borders, inner], to be specific, [bottom, right, up, left, inner]
        % this order is determined by bezier_B
        K  = multi_assemble_k(dofid, material, ke0, ke1, iK_mi, jK_mi);        
        Me = build_M_simply(bdofs_num, K, bezier_B);

%         % [ndofs, v+p]. 
        Re = [bezier_B; Me]; % NOTE: current order of ndofs is as of [dofid]       
        R1 = Re(dofid(order1), :);
        
        % ke = Re' * K * Re;
        % sK(:, ele) = ke(:);
        R2 = reshape(R1, 8, []); % [8, nele*v]
        
        for gp = 1:4
            dNe_g = dNh_ele{gp};
            dFe = dNe_g * R2; % [4, nele*v]
            dF1 = reshape(dFe, 4, eleNum_mi, vdofs_num);
            dFe = permute(dF1, [1,3,2]);
            dFdx{gp, ele} = dFe;
        end
        regulated_matrix{ele} = Re;    
        
%         %% surface of NH
%         R3 = Re(dofid, :);
%         NHe = NH_alpha(microx, edofMat_mi, vdofs_num/2, Nh, R3);
%         NH{ele, 1} = NHe;
    end   

%     plot_NH_bezier(1, 1, lx, microx, nelx, nely, 2, NH);

%     plot_NH(1, 1, lx, microx, nelx, nely, num_b, NH);
    
end