%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-12-16 14:45:35
 % @LastEditTime: 2021-01-10 15:11:37
%--------------------------
function K = update_stiffness(dNh, edofMat_mi, edofMat_ma, n, global_u, global_l, dx, iK, jK)
    eleNum_mi = size(edofMat_mi, 1);
    eleNum_ma = size(edofMat_ma, 1);
    I = eye(eleNum_mi);
    
    jac = dx^2/4;    
    sK = zeros(size(edofMat_ma,2)^2, eleNum_ma);
    
    tmp_edof = edofMat_mi';

    for ele_ma = 1:eleNum_ma
        u = global_u{ele_ma}(:);
        l = global_l{ele_ma}(:);
        
        ne = n{ele_ma}; % [2v, 2n]
        n1 = ne(:, tmp_edof(:))'; % [8m, 2v]

        ke = 0;
        for gp = 1:4
            dNh_g = dNh{gp};        % [4,  8]
            dNh1  = kron(I, dNh_g); % [4m, 8m]
            dNH   = dNh1 * n1;      % [4m, 2v]

            % reshape dNH to [4,2v,m]
            dN1 = reshape(dNH, 4,eleNum_mi,[]); % [4,m,2v]
            dN  = permute(dN1, [1,3,2]);        % [4,2v,m]

            dP1 = dPdF_fast_linear(u, l, dN);    % 4 * l * m
            dPT = permute(dP1, [2,1,3]);         % l * 4 * m 
            dP2 = reshape(dPT, [], 4*eleNum_mi);      % l * 4m           
                        
            dF1 = permute(dN, [2,1,3]);     % l * 4 * m
            dF1 = dF1(:, [1,3,2,4], :); % NOTE:
            dF2 = reshape(dF1, [], 4*eleNum_mi)'; % 4m * l

            ke_g = dP2 * dF2;                % lenout * lenout            
            ke = ke + ke_g * jac;
        end
        sK(:, ele_ma) = ke(:);
    end
    K = sparse(iK, jK, sK);
    K = (K + K')/2;
end