%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-10-29 15:04:18
 % @LastEditTime: 2021-01-12 15:09:02
%--------------------------
function [deform1, converged] = line_search(rhs, du, fext, deform0, ref_nodes, ...
        dN_c, gradN, dF, edofMat, fixeddofs, lx, u, l, opt_rotation)
    maxloop = 10; loop = 0;
    nodeNum = size(deform0, 1);
    deform0 = reshape(deform0', [], 1);
    ref0    = reshape(ref_nodes', [], 1); 

    step = 1;
    approxE0 = -sum(dot(rhs, deform0 - ref0));
    norm0 = max(max(abs(rhs)));
    found_step = 0;

    while loop < maxloop
        loop = loop+1;

        deform1 = step*du + deform0;
        deform2 = reshape(deform1,2,[])';
        U = reshape(deform2' - ref_nodes', [], 1);
        if opt_rotation            
            R = update_corotation_matrix(edofMat, U, dN_c);
            [F, dF] = corotational_F(R, edofMat, ref_nodes, deform2, gradN);
        else
            F = deform_grad_1scale(edofMat, U, gradN); % {4}: 4*eleNum
        end
        % F = deform_gradient(eles, reshape(deform1,2,[])', ref_nodes, gradN); % {4*eleNum}: 2*2

        inverted = check_inversion(F);
        step = step*0.5;
        if (inverted)
            continue; 
            % regard as not converged, 
            % it is advised that the deformation gradient F be temporarily replaced by 
            % the nearest physically plausible value ~F
        end

        fint = elastic_force(F, dF, nodeNum, edofMat, lx, u, l);
        rhs = fint + fext;
        rhs(fixeddofs) = 0;

        approxE1 = -sum(dot(rhs, deform1 - ref0));
        norm1 = max(max(abs(rhs)));
        if( (norm1<norm0) || (approxE1<approxE0) )
            found_step=1;
            break;
        end
    end

    converged = 0; 
    if(~found_step || (abs(norm1-norm0)<1e-5))
        converged = 1;
    end

    deform1(fixeddofs) = deform0(fixeddofs);
    deform1 = reshape(deform1, 2, [])';
end

function inverted = check_inversion(F)
    
    inverted = false;
    eps = 1e-5;
    
    for gp = 1:4
        F1 = F{gp, 1}; % 4 * eleNum
        Fdet = F1(1, :).*F1(4, :) - F1(2, :).*F1(3, :);
        if (sum(Fdet <= eps) > 0)
            inverted = true;
                return;
        end
    end
end