%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-02-22 14:18:49
 % @LastEditTime: 2021-02-22 16:09:07
%--------------------------
% I = eye(3);
% dP = u * (dF + dF') + l * trace(dF) * I;
function dP = dPdF_fast_linear(global_u, global_l, dF)
    
    [partial_num, dofs_num, m] = size(dF); 
    dFT = dF([1:3:9,2:3:9,3:3:9], :, :);     

    %% dP1 = u * (dF + dF')
    mu(1,1,:) = global_u(:);
    dP1 = (dF + dFT) .* mu; % [9,24,m]

    %% dP3 = l * trace(dF) * I;
    dF_expand  = reshape(dF, partial_num, []); % [9,24,m] -> [9,24*m]    
    
    lambda = global_l(:)';               % 1 * m
    trF = reshape(dF_expand(1,:)+dF_expand(5,:)+dF_expand(9,:), dofs_num, []); % 24*m
    lm_tr(1,:,:) = lambda .* trF; % 1*24*m
    
    I  = eye(3);
    I1 = reshape(repmat(I(:), dofs_num, m), size(dF)); % expand eye(3) to 9*24*m
    dP2 = lm_tr .* I1;

    dP = dP1 + dP2;
end