%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-02-22 13:20:56
 % @LastEditTime: 2021-08-24 16:37:34
%--------------------------
function K = multi_stiffness_linear(dFdx, iK, jK, dx, global_u, global_l, lenout, D1, rho_cell)

    eleNum = size(dFdx, 2);    
    jac = dx^2/4;

    sK = zeros(lenout^2, eleNum);

    for ele_ma = 1:eleNum      
%         tic
        u = global_u{ele_ma};
        l = global_l{ele_ma};

        rho_e = rho_cell{ele_ma}(:); % m,1
        tmp = speye(length(rho_e)) .* rho_e; % m,m
        D = kron(tmp, D1);

        % element stiffness
        ke = 0; 
        for gp = 1:4
            dFe = dFdx{gp, ele_ma}; % 4 * lenout * eleNum_mi           

            % dP1 = dPdF_fast_linear(u, l, dFe);   % 4 * l * m
            % dPT = permute(dP1, [2,1,3]);             % l * 4 * m 
            % dP2 = reshape(dPT, lenout, []);          % l * 4m
                        
            dF1 = permute(dFe, [2,1,3]);     % l * 4 * eleNum_mi 
            dF2 = reshape(dF1, lenout, [])'; % 4m * l
                        
            ke_g = dF2' * D * dF2;                % lenout * lenout            

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