%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-05-18 19:11:25
 % @LastEditTime: 2021-05-18 21:54:48
%--------------------------
% fast version
% dF.size() = 9 * 24 * m
function dP = dPdF_fast_nonlinear(global_u, global_l, F, dF)   

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

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

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

function c = matrix_3d_multiply_3rank(a, b)
    % [ a1*b1 + a4*b2 + a7*b3, a1*b4 + a4*b5 + a7*b6, a1*b7 + a4*b8 + a7*b9]
    % [ a2*b1 + a5*b2 + a8*b3, a2*b4 + a5*b5 + a8*b6, a2*b7 + a5*b8 + a8*b9]
    % [ a3*b1 + a6*b2 + a9*b3, a3*b4 + a6*b5 + a9*b6, a3*b7 + a6*b8 + a9*b9]
    c = [a(1,:) .* b(1,:) + a(4,:) .* b(2,:) + a(7, :) .* b(3, :);
        a(2,:) .* b(1,:) + a(5,:) .* b(2,:) + a(8, :) .* b(3, :);
        a(3,:) .* b(1,:) + a(6,:) .* b(2,:) + a(9, :) .* b(3, :);
        a(1,:) .* b(4,:) + a(4,:) .* b(5,:) + a(7, :) .* b(6, :);
        a(2,:) .* b(4,:) + a(5,:) .* b(5,:) + a(8, :) .* b(6, :);
        a(3,:) .* b(4,:) + a(6,:) .* b(5,:) + a(9, :) .* b(6, :);
        a(1,:) .* b(7,:) + a(4,:) .* b(8,:) + a(7, :) .* b(9, :);
        a(2,:) .* b(7,:) + a(5,:) .* b(8,:) + a(8, :) .* b(9, :);
        a(3,:) .* b(7,:) + a(6,:) .* b(8,:) + a(9, :) .* b(9, :)];    
end