%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-11-28 21:26:00
 % @LastEditTime: 2021-02-28 18:45:51
%--------------------------
function [U1, converged] = line_search_linear(rhs, du, U0, fext, dN, edofMat, fixeddofs, lx, u, l)
    
    maxloop = 10; loop = 0;
    alldofs = size(fext, 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;
        F = deform_grad(dN, edofMat, U1);

        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, dN, alldofs, 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);
end

function inverted = check_inversion(F)

    inverted = false;
    eps = 1e-5;

    for gp = 1:8
        F1 = F{gp, 1}; % 9 * eleNum
        Fdet = det_3rank(F1);
        if (sum(Fdet <= eps) > 0)
            inverted = true;
            return;
        end
    end
end