% To check and understand Main_Xarxa used for first submission
% 21/09/2021
clear all
close all
addpath Functionals:Interpolation_Geometry:Reparametrization:Examples_static_NEW:Matrices
global phi Invariant_aux C_1 dt ini_Val Material Constraints V CC  C_psi 


input_file ='Temp_Peak_Height_Vol.mat';
load(input_file)


MaxFunEvals=3000;
MaxIter=50000;
tol_X=1e-3;

% velocity of growth of dendritic network
v0_dendrite = 9e-3;
F0          = 0.5; %0.308
decay_f = 5;
% Specifies whether the dendritic network is in contact with the membrane
contact = 0;
% Change in magnitude of height (just to initialize here)
Ht_Chng = 0;

time = 0;
dt = 1;

Proj_area_0 = pi * x(end-1)^2;
Areal_Strain = 1;
max_psi_hist = [];
max_phi_hist = [];
max_V_c_hist = [];
tip_hist       = [];
tip_force_hist = [];

Material.R_ini = x(end-1);



% options = optimoptions('fmincon','Algorithm','interior-point' ,'Display','final', 'CheckGradient', false, ...
%     'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,'MaxFunctionEvaluations',MaxFunEvals);
% 
options = optimoptions(options,'Algorithm','interior-point',"EnableFeasibilityMode",true,...
        'TolX',1e-7,'Display','final');



C_phi = zeros(size(C_phi));
psi = zeros(size(psi));
%%
while time(end) < t_f & i < i_max & k < 100000
    i = i+1;
    k = k+1;
    fprintf('Instant %10.5f towards %10.5f\n', time(end)+dt, t_f)
        
    %
    [C,C1,C2] = interpolate(Invariant,x);
    Aux = auxiliary_calculations(Invariant,ini_Val.p,Constraints.C0,C,C1,C2);
    %     Constraints.S0_L = localArea_functional(Aux);
    %     Constraints.S0_L = Constraints.S0_L * Constraints.S0 / sum(Constraints.S0_L);
    
    
    max(abs(Aux.dau))
    
    if any(abs(Aux.dau) > max_dau)
        disp('**********************')
        disp('We reparametrize')
        [C_fit,phi_fit] = original_shape_to_fit(x, C_phi); % at points uniformly spaced in arc-length
        %[C_fit] = original_shape_to_fit(x); % at points uniformly spaced in arc-length
        ce=ini_Val.ce; cc=ini_Val.cc;
        %creparametrize(x,ce,cc);
        x_new = reparametrize(C_fit,x);
        x=x_new;
        % Update protein coverage
        C_phi_new = reparametrize_func(phi_fit,C_phi);
        C_phi=C_phi_new;
        [phi, d_phi, dd_phi] = interpolate_func_2(Invariant,C_phi);

        
        % Update the local area constraints
        [C,C1,C2] = interpolate(Invariant,x);
        Aux = auxiliary_calculations(Invariant,ini_Val.p,Constraints.C0,C,C1,C2);
        Constraints.S0_L = localArea_functional(Aux);
        Constraints.S0_L = Constraints.S0_L * Constraints.S0 / sum(Constraints.S0_L);
        ct_rep=ct_rep+1;
        if(ct_rep==2)
            max_dau=max_dau*1.15;
            if(max_dau>200)
                max_dau = 200;
            end
            ct_rep=0;
        end
    else
        ct_rep=0;
    end
    
    phi_1 = phi;
    C_phi_1 = C_phi;
    if i>0
        C_phi_k = C_phi;
        [phi_k, d_phi_k, dd_phi_k] = interpolate_func_2(Invariant,C_phi_k);
        R = Residual_mod_Ener_Aug(Aux,C,C1,C2,phi_k,d_phi_k,phi_1,dd_phi_k);
        norm(R);
        M=matrix_Mass(Aux,C,C1,C2);
        A=matrix_AP2(Aux,C,C1,C2);

        ctr=1;
        while (norm(R) > 1e-10 && ctr <10)
            K=matrix_K_mod_Ener_Aug(Aux,C,C1,C2,phi_k);
            B=matrix_B_mod_Ener_Aug(Aux,C,C1,C2,phi_k,d_phi_k);
            JD=matrix_JD(Aux,C,C1,C2,phi_k,d_phi_k,dd_phi_k);
            J = (M + A + K + B +JD);
            J(:,end) = 0;
            J(end,:) = 0;
            J(end,end)=1;
            R(end)=0;
                     
            C_phi_new = C_phi_k - J'\R;
            C_phi_k= C_phi_new;
            [phi_k, d_phi_k, dd_phi_k] = interpolate_func_2(Invariant,C_phi_k);
            R = Residual_mod_Ener_Aug(Aux,C,C1,C2,phi_k,d_phi_k,phi_1,dd_phi_k);
%            R = Residual_IncStiff(Aux,C,C1,C2,phi_k,d_phi_k,phi_1,dd_phi_k);
            ctr=ctr+1;
        end
        [phi_k, d_phi_k, dd_phi_k] = interpolate_func_2(Invariant,C_phi_k);
        C_phi = C_phi_k;
        phi = phi_k;
        d_phi = d_phi_k;
        dd_phi = dd_phi_k;
    end
    
    
    if 1==1
        % Map phi to the flat cortex
        [N ,EDGES ,BIN ] = histcounts(C(:,1),ini_Val.UU*Material.R0);
        [N_,EDGES_,BIN_] = histcounts(Invariant.gaussp,ini_Val.UU);
        phi_flat_ = zeros(size(Invariant.gaussp));
        for kk = 1:length(ini_Val.UU)-1
            if sum(BIN==kk)>0
                phi_flat(kk) = max(phi(BIN==kk));
            else
                phi_flat(kk)=phi_flat(kk-1);
            end
            phi_flat_(BIN_==kk) = phi_flat(kk);
        end
        %figure(77)
        %hold off
        %plot(0.5*(ini_Val.UU(1:end-1)+ini_Val.UU(2:end)), phi_flat, 'ko')
        %hold on
        %plot(Invariant.gaussp, phi_flat_, 'rx')
        
        
        psi_1 = psi;
        C_psi_1 = C_psi;
        
        M_psi=matrix_Mass_psi;
        K_psi=matrix_Stiffness_psi;
        Rhs_psi = vector_rhs_psi(psi_1,phi);
        C_psi = (M_psi + K_psi)\Rhs_psi;
        [psi, d_psi] = interpolate_func_2(Invariant,C_psi);

%         figure(78)
%         plot(Invariant.gaussp, psi)
%         xlabel('r')
%         ylabel('psi')
        Cortex_Flat_MAB(psi);
    end
    
    
    x_1 = x;
    [C_1,C1_1,C2_1] = interpolate(Invariant,x_1);
    [x, fval, exitflag, output, lambda] = fmincon(@Incremental_Mech_Rayleighian_S,x,[],[],Aeq,beq,[],[],@Nonlin_constr,options);
    
    
    if i>0
        if  max(abs(x-x_1)) > 0.5
            disp('*****************')
            disp('We reduce the timestep')
            dt = dt/2;
            if(contact)
                Ht_Chng  = Ht_Chng/2;
                beq(4)   = beq(4) - Ht_Chng;
            end
            pause(0.2)
            x = x_1;
            continue
        end
        if max(abs(x-x_1)) < 5e-2
            dt = dt * 1.5;
            if dt >= dt_max
                dt = dt_max;
                
            else
                disp('*****************')
                disp('We increase the timestep')
            end
        end
        
    end
    
    % To change the height depending on the force

    

    Material.eta_END = Material.eta_s * pi * x(end-1)^2 * log(100);

    [C,C1,C2] = interpolate(Invariant,x);
    Aux = auxiliary_calculations(Invariant,ini_Val.p,Constraints.C0,C,C1,C2,phi);
    
    Invariant_mem_to_cortex=precal_generic(C(:,1)/Material.R0);

    %[psi_m]=interpolate_func_2(Invariant_mem_to_cortex,C_psi);
    %[tD_c,dD_dCP_c] = Incremental_cortex_diss(Aux,C,C_1,C1);
     
    if(contact)
        disp(lambda.eqlin(4))
        Ht_Chng = dt*v0_dendrite*(exp(decay_f*lambda.eqlin(4)) - exp(-decay_f*F0));
        beq(4)  = beq(4) + Ht_Chng;
    end
    
    
    time = [time time(end)+dt];
    tip_hist = [tip_hist x(2)];
    
    if(contact)
        tip_force_hist = [tip_force_hist lambda.eqlin(4)];
    else
        tip_force_hist = [tip_force_hist 0];
    end
    
    C_psi_hist = [C_psi_hist C_psi];
    [phi_s] = interpolate_func_2(Invariant_sample,C_phi);
    phi_s_hist = [phi_s_hist phi_s];
    [C_s,C1_s,C2_s] = interpolate(Invariant_sample,x);
    C_s_hist = [C_s_hist C_s];
    V_hist = [V_hist V];



    area_hist = [area_hist Area_functional(Aux)];
    max_psi_hist  = [max_psi_hist  max(psi)];
    max_phi_hist  = [max_phi_hist  max(phi_s)/Material.phi_bar];

    
    
    if (Constraints.iSurf_loc ==1) %&& count==0
        tension = interpolate_func(Invariant,lambda.eqnonlin(2:end),Invariant.lower);
        
        
    else
        tension = lambda.eqnonlin*ones(size(Invariant.gaussp));
        %count=1;
    end
    sigma = [sigma tension];
    
    %plot respect to archlength
    D= (Aux.au.* Invariant.gweight);
    S=zeros(size(D));
    for i=1:size(D)
        S(i)= sum(D(1:i));
    end
    %%%%%%%%%%%%%%%%%%%%%
    figure(3)
    hold off
    plot(C_s( :,1),C_s( :,2))
    hold on
    plot(x(1:2:end),x(2:2:end),'ko')
    axis equal
    %Cortex velocity on membrane
    [v_c_m]=interpolate_func_2(Invariant_mem_to_cortex,V);
    [v_c_c]=interpolate_func_2(Invariant,V);
    quiver(Invariant.gaussp*Material.R0,Invariant.gaussp*Material.R0*0,v_c_c',v_c_c'*0)
    tmp_z  = 0.5*(1-tanh(Material.tmp_s*(C(:,2)-Material.h_cut_off)));
    v_c_m = v_c_m .* tmp_z;
    quiver(C(:,1),C(:,2),v_c_m,v_c_m*0)
    %quiver(x(1:2:end), x(2:2:end), dD_dCP_c(1,:)', dD_dCP_c(2,:)') 
    
    
    max_V_c_hist  = [max_V_c_hist  max(v_c_c)];
   
    
     if(v0_dendrite*(1-exp(-decay_f*F0))*time(end)>x(2)  && contact==0)
        dt          = dt/5; 
        Aeq         = [Aeq; zeros(1,length(Aeq))];
        Aeq(4,2)    = 1;
        Ht_Chng     = dt*v0_dendrite;
        beq(4)      = x(2) + Ht_Chng;
        contact     = 1;
        disp('We establish contact between growing dendritic network and membrane')
    end
    
    
    Proj_area_i = pi*x(end-1)^2;
    %Areal_Strain = [Areal_Strain 100.*(Proj_area_i - Proj_area_0)./Proj_area_0];
    Areal_Strain = [Areal_Strain (pi*Material.R0^2 - Proj_area_i )/(pi*Material.R0^2 - Proj_area_0 )];
    figure (7)
    title('Normalized excess area vs Time')
    plot(time,Areal_Strain,'--ro','Markersize',4,'LineWidth',2)
    xlabel('Time','FontSize',20,'FontWeight','bold','Color','k')
    ylabel('Normalized excess area','FontSize',20,'FontWeight','bold','Color','k')
    pause(0.1)
    
%     figure (8)
%     title('Density vs position')
%     plot(C(:,1),phi,'--bo','Markersize',4,'LineWidth',2)
%     %plot(Invariant.gaussp,phi,'--bo','Markersize',4,'LineWidth',2)
%     xlabel('R','FontSize',20,'FontWeight','bold','Color','k')
%     ylabel('Density','FontSize',20,'FontWeight','bold','Color','k')
%     pause(0.1)
    
end
%C_phi_hist = [C_phi_hist_1 C_phi_hist];

figure(990)
plot(time(2:end),max_phi_hist)
xlabel('time')
ylabel('maxphi')
figure(991)
plot(time(2:end),max_psi_hist)
xlabel('time')
ylabel('maxpsi')
figure(992)
plot(time(2:end),max_V_c_hist)
xlabel('time')
ylabel('maxVc')

return
%% 
output_file = 'Examples_Dynamic/Bud_Eq_F0_0_12_Dendrite_V0_9e3_F0_0_22_Exp_D_25_No_Cortex_Flow.mat';

osave = 1;
if osave == 1
    save(output_file);
end



