%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-11-28 21:26:00
 % @LastEditTime: 2021-01-12 15:10:03
%--------------------------
function [U1, converged] = line_search_linear(rhs, du, fext, U0, gradN, dF, edofMat, fixeddofs, lx, u, l)
    maxloop = 10; loop = 0;
    nodeNum = size(fext, 1) / 2;
%     deform0 = reshape(deform0', [], 1);
%     ref0    = reshape(ref_nodes', [], 1); 

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

    while loop < maxloop
        loop = loop+1;

        U1 = step*du + U0;
%         deform2 = reshape(deform1,2,[])';
        F = deform_grad_1scale(edofMat, U1, gradN);
        
%         if opt_rotation
%             U = reshape(deform2' - ref_nodes', [], 1);
%             R = update_corotation_matrix(edofMat, U, dN_c, deform2);
%             [F, dF] = corotational_F(R, edofMat, ref_nodes, deform2, gradN);
%         else
%             F = deform_gradient(eles, deform2, ref_nodes, 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_linear(F, dF, nodeNum, edofMat, lx, u, l);
        rhs = fint + fext;
        rhs(fixeddofs) = 0;

        approxE1 = -sum(dot(rhs, U1));
        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

    U1(fixeddofs) = U0(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