%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-02-28 16:50:32
 % @LastEditTime: 2021-03-02 09:56:36
%--------------------------
function K = stiffness_2scale_linear(dFdx, iK, jK, dx, global_u, global_l, vdofs_num)

    nele = size(dFdx, 2);    
    jac = dx^3 / 8;

    sK = zeros(vdofs_num^2, nele);

    parfor ele_ma = 1:nele        
        u = global_u{ele_ma};
        l = global_l{ele_ma};

        % element stiffness
%         tic
        ke = 0; 
        for gp = 1:8
            dFe = dFdx{gp, ele_ma}; % 9 * v * eleNum_mi           

            dP1 = dPdF_fast_linear(u, l, dFe);   % 9, v, m
            dPT = permute(dP1, [2,1,3]);         % v, 9, m
            dP2 = reshape(dPT, vdofs_num, []);   % v, 9m
                        
            dF1 = permute(dFe, [2,1,3]);         % v, 9m
            dF2 = reshape(dF1, vdofs_num, [])';  % 9m * v
                        
            ke_g = dP2 * dF2;                % v * v            
            % ke_g = sum(mtimesx(dPT, dFe), 3); % 24*24*m
            ke = ke + ke_g * jac;
        end
        sK(:, ele_ma) = ke(:);
%         toc
    end
    K = sparse(iK, jK, sK(:));
    K = (K + K')/2;
end