%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-10-29 15:04:18
 % @LastEditTime: 2021-01-11 22:00:36
%--------------------------
% fast version
% dF.size() = 4 * l * m
function dP = dPdF_fast(global_u, global_l, F, dF)    

    % expand & reorder -> (4, lenout*eleNum_mi)
    lenout = size(dF,2);
    F1  = repmat(F, lenout, 1);
    F1 = permute(F1, [1,3,2]);
    F1  = reshape(F1(:), 4, []);
    dF_expand = reshape(dF, 4, []);

    dFT = dF_expand([1,3,2,4], :); % NOTE: change 2 rows to get dF^T

    Fdet = F1(1, :).*F1(4, :) - F1(2, :).*F1(3, :);
    Finv = F1([4,2,3,1], :) .* [1;-1;-1; 1] ./ Fdet;    

    %% dP1 = u * dF
    mu(1,1,:) = global_u(:);
    dP1 = dF .* mu; % 4 * l * m
    
    %% dP2 = (mu - lambda * JJ) * Finv' * dF' * Finv';
    JJ = reshape(log(Fdet), lenout, []); % l * m
    lambda = global_l(:)';               % 1 * m
    lm_J(1,:,:) = lambda .* JJ;          % 1 * l * m
        
    FinvT = Finv([1,3,2,4], :);
    tmp = matrix_3d_multiply(matrix_3d_multiply(FinvT, dFT), FinvT);    
    tmp = reshape(tmp, size(dF));
    
    dP2 = mu .* tmp - lm_J .* tmp;
    
    %% dP3 = lambda * trace(F\dFe) * Finv'
    tmp = matrix_3d_multiply(Finv, dF_expand);
    trF = reshape(tmp(1,:)+tmp(4,:), lenout, []); % l*m
    lm_tr(1,:,:) = lambda .* trF; % 1*l*m
    
    dP3 = lm_tr .* reshape(FinvT, size(dF));

    dP = dP1 + dP2 + dP3; % 4 * l * m     
end