%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-11-28 21:07:50
 % @LastEditTime: 2020-11-28 21:49:42
%--------------------------
function f = elastic_force_linear(F, dF, nodeNum, edofMat, lx, u, l)

    f = zeros(nodeNum*2, 1);
    jac = lx^2/4;
    eleNum = size(edofMat,1);
    
    for gp = 1:4
        dFe = dF{gp};                      % 4 * 8 * m
        Fe  = F{gp, 1};                    % 4 * m
        P2(1, :, :)  = PK1_fast_linear(u, l, Fe); % 4 * m
        
        fe = squeeze(mtimesx(P2, dFe)) * jac; % 8*m

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

% I = eye(2,2);
% PP = u*(F + F'-2*I) + l * trace(F-I) * I;
function P = PK1_fast_linear(u, l, F)
    m = size(F,2);
    
    FT  = F([1,3,2,4], :);     % 4*m
    trF = F(1,:) + F(4,:) - 2; % 1*m

    I = repmat([1;0;0;1], 1, m); % 4*m

    P = u(:)' .* (F + FT - 2*I) + l(:)' .* trF .* I;
end