%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-11-10 21:04:46
 % @LastEditors: Please set LastEditors
 % @LastEditTime: 2020-11-11 14:16:26
%--------------------------
function f = elastic_force_homog(lenout, nodeNum, dx, F, dFdx, edofMat_ma, u0, l0)
    eleNum_ma = size(dFdx, 2);
    f = zeros(nodeNum*2, 1);
    jac = dx^2/4;

    for ele_ma = 1:eleNum_ma
        
        % as micro outer order: 4_coners, bottom, right, up, left
        fe = 0;
        for gp = 1:4
            Fe = F{gp, ele_ma};      % 4 * eleNum_mi
            P2 = PK1_homog(u0, l0, Fe); % 4 * eleNum_mi

            dFe = dFdx{gp, ele_ma};          % 4 * lenout * eleNum_mi
            dF1 = permute(dFe, [2,1,3]);     % l * 4 * eleNum_mi 
            dF2 = reshape(dF1, lenout, [])'; % 4m * l

            fg  = P2(:)' * dF2;  % 1 * lenout

            fe = fe + fg * jac;         
        end

        % every macro ele: 4corners -> down -> right -> up -> left
        edof = edofMat_ma(ele_ma, :); 
        f(edof) = f(edof) - fe(:); % NOTE: MINUS!
    end   
end

function P = PK1_homog(u, l, F)
    Fdet = F(1, :).*F(4, :) - F(2, :).*F(3, :);
    JJ   = log(Fdet);
    Finv = F([4,2,3,1], :) .* [1;-1;-1; 1] ./ Fdet;
    FinvT = Finv([1,3,2,4], :);

    P = u * (F - FinvT) + l * JJ .* FinvT;
end