%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-10-29 15:04:18
 % @LastEditTime: 2021-02-22 15:40:56
%--------------------------
function K = multi_global_stiffness(F, dFdx, iK, jK, dx, global_u, global_l, lenout)

    eleNum = size(F, 2);
    sK = zeros(lenout^2, eleNum);
    jac = dx^2/4;

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

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

            % NOTE: dP = function(dF), and dF of every dofs is [dNdx, dNdy; 0, 0] or [0, 0; dNdx, dNdy]
            dP1 = dPdF_fast(u, l, Fe, 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

            % NOTE: original dP * dF is [2,2]*[2,4], i.e. [dp1, dp2;   [dNdx
            %                                              dp3, dp4] *  dNdy]
            % here is [1,4]*[4,8], i.e. [dp1, dp3, dp2, dp4] * [dNdx; 0; dNdy; 0; ...]
            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