%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-02-22 16:27:54
 % @LastEditTime: 2021-02-22 19:40:29
%--------------------------
function f = elastic_force_linear(F, dNe, alldofs, edofMat, lx, u, l)

    f = zeros(alldofs, 1);
    jac = lx^3/8;
    eleNum = size(edofMat,1);
    
    for gp = 1:8
        dNe_g = dNe{gp};              % 9,24
        dFe = repmat(dNe_g, [1,1,eleNum]); % 9,24,m
        Fe  = F{gp, 1};               % 9 * m
        P2(1, :, :)  = PK1_fast_linear(u, l, Fe); % 1 * 9 * m
        
        fe = squeeze(mtimesx(P2, dFe)) * jac; % 24*m

        for ele = 1:eleNum   
            edof = edofMat(ele, :);
            f(edof) = f(edof) - fe(:, ele); % NOTE: MINUS!
        end        
    end 
end