%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-11-28 21:26:00
 % @LastEditTime: 2021-05-18 20:52:06
%--------------------------
function [U1, converged] = line_search_2scale_linear(rhs, du, U0, fext, edofMat, fixeddofs, lx, u, l, ...
    sort_u, dFdx, dNe, nele_mi, regulated_matrix)
    
    maxloop = 10; loop = 0;
    alldofs = size(fext, 1);
    vdofs_num = size(edofMat, 2);

    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_2scale(sort_u, edofMat, dNe, nele_mi, U1, regulated_matrix);

        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_2scale_linear(vdofs_num, alldofs, lx, F, dFdx, edofMat, 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