clear all

%%

% Charger les données d'un fichier .mat
filename = 'Test8.mat';
load(filename)

% Obtenir le nombre de noeuds et d'éléments dans le modèle
NNode = length(Model.Nodes.XNOD);
NElem = length(Model.Elements.ELEMNOA);

% Définir les degrés de liberté bloqués
listDDLbloque = [1 2 3 3*NNode-2 3*NNode-1 3*NNode];
lim_graph = [4 9];
% Définir une série de constantes de raideur pour les ressorts
K_ressort = [2 4 6 8 10 15 20 25 40 60 80 100 1000 1e4];
kappa = zeros(size(K_ressort));
% Initialiser les tableaux pour stocker les résultats
vect_L1 = zeros(length(K_ressort),1);
vect_L_FEM = zeros(length(K_ressort),1);
vect_sum = zeros(length(K_ressort),1);
vect_H_sigma = zeros(length(K_ressort),1);
e = zeros(length(K_ressort),1);
vect_L1P = zeros(length(K_ressort),1);
vect_sumP = zeros(length(K_ressort),1);
Rapport = zeros(length(K_ressort),1);

% Cas infiniment rigide
Model.Springs(3).Value = 1e18;

% Calcul des modes propres et de la matrice de rigidité dans le cas
% infiniment rigide
[K0, Ksigma, phi0, L_FEM, Re, Model, ddlLibres] = calcul_FEM(Model,listDDLbloque);
K00 = K0(ddlLibres,ddlLibres);
Ksigma00 = Ksigma;
L0 = L_FEM;
%% Calcul de Lambda1 avec la première méthode

vect_L1_elem = zeros(NElem,length(K_ressort));

for i = 1:length(K_ressort)
  
    % Calcul de Lambda1
    DH0_elem = 0;
    DHsigma_elem = 0;
    Hsigma_elem = 0;
    DH0 = 0;
    DHsigma = 0;
    Hsigma = 0;

    % Parcourir tous les éléments pour calculer Lambda1
    for j = 1:NElem
        %Données
        I_typ = Model.Elements.ELEMGEO(j);
        I = Model.Geometries(I_typ).I;
        E = Model.Materials(1).Young;
        EI = E*I;
        
        NOA = Model.Elements.ELEMNOA(j);
        NOB = Model.Elements.ELEMNOB(j);
        Le = Model.Elements.ELEMLEN(j);
        P = Model.Elements.N_STAB(j);
        phi_elem(1:3) = Re(1:3,1:3,j)'*phi0(3*NOA-2:3*NOA);
        phi_elem(4:6) = Re(4:6,4:6,j)'*phi0(3*NOB-2:3*NOB);
        phi1 = phi_elem(3);
        phi2 = phi_elem(6);
        

        % Calculer les coefficients alpha pour les ressorts
        if Model.Elements.ELEMSP1(j) == 3
            Model.Springs(3).Value = K_ressort(i);
            alpha1 = EI/(K_ressort(i)*Le);
        end
        if Model.Elements.ELEMSP2(j) == 3
            Model.Springs(3).Value = K_ressort(i);
            alpha2 =  EI/(K_ressort(i)*Le);
        end
        if Model.Elements.ELEMSP1(j) ~= 3
            alpha1 = 0;
        end
        if Model.Elements.ELEMSP2(j) ~= 3
            alpha2 = 0;
        end

            psi = (phi_elem(5) - phi_elem(2))/Le;

            %Calcul de DH0
            DH0_elem = -4*EI/Le*(alpha1*(2*phi1 + phi2 - 3*psi)^2 + alpha2*(phi1 + 2*phi2 - 3*psi)^2);
            DH0 = DH0 - 4*EI/Le*(alpha1*(2*phi1 + phi2 - 3*psi)^2 + alpha2*(phi1 + 2*phi2 - 3*psi)^2);

            %Calcul de DHsigma
            DHsigma_elem = - 2*P*Le/15*(alpha1*(2*phi1 + phi2 - 3*psi) * (4*phi1 - phi2 - 3*psi) ...
                    + alpha2*(phi1 + 2*phi2 - 3*psi)*(4*phi2 - phi1 - 3*psi));
            DHsigma = DHsigma - 2*P*Le/15*(alpha1*(2*phi1 + phi2 - 3*psi) * (4*phi1 - phi2 - 3*psi) ...
                    + alpha2*(phi1 + 2*phi2 - 3*psi)*(4*phi2 - phi1 - 3*psi));

            %Calcul de Hsigma
            Hsigma_elem = P*Le/15 * (2*phi1^2 - phi1*phi2 + 2*phi2^2 - 3*(phi1+phi2)*psi + 18*psi^2);
            Hsigma = Hsigma + P*Le/15 * (2*phi1^2 - phi1*phi2 + 2*phi2^2 - 3*(phi1+phi2)*psi + 18*psi^2);
            
            vect_L1_elem(j,i) = -(DH0_elem + L0*DHsigma_elem)/Hsigma_elem;
    end

        
        L1 = -(DH0 + L0*DHsigma)/(Hsigma);

        Lambda_sum = L0 - L1;
        
        % Calcul du L_FEM
        [~,~, ~, L_FEM, ~, Model,ddlLibres] = calcul_FEM(Model,listDDLbloque);
        

% On stock les valeurs importantes

vect_L1(i) = L1;
vect_rapport(i) = L1/L0;
vect_L_FEM(i) = L_FEM;
vect_sum(i) = Lambda_sum;
e(i) = (L_FEM - Lambda_sum)/L_FEM * 100;
kO(:,:,i) = K0;
kappa(i) = K_ressort(i)*Le/(Model.Geometries(end).I);


end


%% Calcul de Lambda1 avec la deuxième méthode
P1 = zeros(length(phi0),1);
Psigma = zeros(length(phi0),1);
phi0_approx = zeros(size(phi0));

F = zeros(size(phi0));

F(13) = 1;
phi0_approx(ddlLibres) = -K00\F(ddlLibres);
phi0_approx = phi0_approx/max(abs([phi0_approx(1:3:end); phi0_approx(2:3:end)]));
%phi0_approx = phi0;

vect_P1 = zeros(6,NElem);
eps_lambda1 = zeros(NElem,1);

for i = 1:length(K_ressort)
    
    Hsigma = 0;

    % Calcul de Lambda1
    P1 = zeros(length(phi0),1);
    Psigma = zeros(length(phi0),1);
    
        % Parcourir tous les éléments pour calculer Lambda1
    for j = 1:NElem
        %Données
        I_typ = Model.Elements.ELEMGEO(j);
        I = Model.Geometries(I_typ).I;
        E = Model.Materials(1).Young;
        EI =E*I;
        
        NOA = Model.Elements.ELEMNOA(j);
        NOB = Model.Elements.ELEMNOB(j);
        Le = Model.Elements.ELEMLEN(j);
        P = Model.Elements.N_STAB(j);
        phi_elem(1:3) = Re(1:3,1:3,j)'*phi0_approx(3*NOA-2:3*NOA);
        phi_elem(4:6) = Re(4:6,4:6,j)'*phi0_approx(3*NOB-2:3*NOB);
        phi1 = phi_elem(3);
        phi2 = phi_elem(6);
        P1_elem = zeros(6,1);
        Psigma_elem = zeros(6,1);

        % Calculer les coefficients alpha pour les ressorts
        if Model.Elements.ELEMSP1(j) == 3
            alpha1 = EI/(K_ressort(i)*Le);
        end
        if Model.Elements.ELEMSP2(j) == 3
            alpha2 = EI/(K_ressort(i)*Le);
        end
        if Model.Elements.ELEMSP1(j) ~= 3
            alpha1 = 0;
        end
        if Model.Elements.ELEMSP2(j) ~= 3
            alpha2 = 0;
        end

        psi = (phi_elem(5) - phi_elem(2))/Le;
        
        
        PHI1 = 2*phi1 + phi2 - 3*psi;
        PHI2 = phi1 + 2*phi2 - 3*psi;
        PHI1_etoile = 3*PHI1 + 2*(phi1 -phi2);
        PHI2_etoile = 3*PHI2 + 2*(phi2 - phi1);
        
        % Compute P1
        P1_elem(1) = 0;
        P1_elem(2) = -12*EI/Le^(2)*(alpha1*PHI1 + alpha2*PHI2);
        P1_elem(3) = -4*EI/Le*(2*alpha1*PHI1 + alpha2*PHI2);
        P1_elem(4) = 0;
        P1_elem(5) = 12*EI/Le^(2)*(alpha1*PHI1 + alpha2*PHI2);
        P1_elem(6) = -4*EI/Le*(alpha1*PHI1 + 2*alpha2*PHI2);
        
        vect_P1(:,j) = P1_elem;
        
    
        P1(3*NOA-2:3*NOA,1) = P1(3*NOA-2:3*NOA,1) + Re(1:3,1:3,j)*P1_elem(1:3);
        P1(3*NOB-2:3*NOB,1) = P1(3*NOB-2:3*NOB,1) + Re(4:6,4:6,j)*P1_elem(4:6);
        
        % Compute Psigma
        Psigma_elem(1) = 0;
        Psigma_elem(2) = -6*P/5*(alpha1*(phi1-psi) + alpha2*(phi2-psi));
        Psigma_elem(3) = -2*P*Le/15*(alpha1*PHI1_etoile - alpha2*(phi1 - phi2));
        Psigma_elem(4) = 0;
        Psigma_elem(5) = 6*P/5*(alpha1*(phi1-psi) + alpha2*(phi2-psi));
        Psigma_elem(6) = -2*P*Le/15*(alpha1*(phi1 - phi2) + alpha2*PHI2_etoile);
        
        Psigma(3*NOA-2:3*NOA,1) = Psigma(3*NOA-2:3*NOA,1) + Re(1:3,1:3,j)*Psigma_elem(1:3);
        Psigma(3*NOB-2:3*NOB,1) = Psigma(3*NOB-2:3*NOB,1) + Re(4:6,4:6,j)*Psigma_elem(4:6);
        
        %Compute Hsigma
        
        Hsigma = Hsigma + P*Le/15*(2*phi1^(2) - phi1*phi2 + 2*phi2^(2) - 3*(phi1 + phi2)*psi + 18*psi^2);
            
        eps_lambda1(j,1) =  P1_elem'*phi_elem';   
            
    end
    
    P1_data = P1;
    Psigma_data = Psigma;

%     for n = 1:length(P1)
%         if mod(n,3) ~= 0
%             P1(n) = 0;
%         end
%     end
    
    P1 = P1(ddlLibres,1);
    Psigma = Psigma(ddlLibres,1);
  

    L1P = (P1')*phi0_approx(ddlLibres)/(Hsigma);
    sumP = L0 + L1P;
    
    vect_L1P(i) = L1P;
    vect_sumP(i) = sumP;
    Rapport(i) = -L1P/L0;
   
end

eps_lambda1 = eps_lambda1/Hsigma;
    

%o = (vect_sumP(7) - vect_sum(7))*100/vect_sum(7)
%% Tracer le graphique

% Calcul de l'erreur en pourcentage entre les deux courbes
err_pct = abs(vect_sum - vect_L_FEM) ./ vect_L_FEM * 100;

% Recherche de l'indice correspondant à l'erreur maximale
[~, idx] = max(err_pct);


figure('Position', [100 100 600 400]);
semilogx(K_ressort, vect_sum, 'b-', 'LineWidth', 1.5, 'MarkerSize', 8);
hold on;
semilogx(K_ressort, vect_sumP, 'k-', 'LineWidth', 1.5, 'MarkerSize', 8);
semilogx(K_ressort, vect_L_FEM, 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 8);
grid on;
legend('$\lambda_{approx,1}$','$\lambda_{approx,2}$', '$\lambda_{FEM}$', 'Interpreter', 'latex', 'Location', 'southeast', 'FontSize', 12);
xlabel('Constante de raideur du ressort (Nm)', 'Interpreter', 'latex', 'FontSize', 12);
ylabel('Multiplicateur critique $\lambda_{cr}$', 'Interpreter', 'latex', 'FontSize', 12);
%title('Relation entre le multiplicateur critique et la constante de raideur du ressort', 'FontSize', 14);
set(gca, 'FontSize', 12);
box on;
ylim(lim_graph);

figure('Position', [700 100 600 400]);
semilogx(K_ressort, Rapport*100, 'k-', 'LineWidth', 1.5, 'MarkerSize', 8);
hold on;
semilogx(K_ressort, vect_rapport*100, 'b-', 'LineWidth', 1.5, 'MarkerSize', 8);
x_intersection = interp1(vect_rapport*100, K_ressort, 5);
x2_intersection = interp1(vect_rapport*100, K_ressort, 1);
grid on;
xlabel('Constante de raideur du ressort (Nm)', 'Interpreter', 'latex', 'FontSize', 12);
ylabel('Rapport $\varepsilon\lambda_{1}/\lambda_{0}$[\%]', 'Interpreter', 'latex', 'FontSize', 12);
ylim([0 20]);

% Ajout de la ligne horizontale en pointillés rouges
yValue = 5;
xRange = [1 x_intersection];
line(xRange, [yValue yValue], 'LineStyle', '--', 'Color', 'r');
line([x_intersection, x_intersection], [0, yValue], 'Color', 'red', 'LineStyle', '--')

% Afficher la valeur de x_intersection sur l'axe des abscisses
text(x_intersection, 6, sprintf('S_{j,ini} = %.2f', x_intersection), 'Color', 'black', 'VerticalAlignment', 'top', 'HorizontalAlignment', 'left')

legend('$\varepsilon\lambda_{1}/\lambda_{0}$(Hypothese)', '$\varepsilon\lambda_{1}/\lambda_{0}$','$\varepsilon\lambda_{1} = 0.05\lambda_{0}$', 'Interpreter', 'latex', 'Location', 'northeast', 'FontSize', 12);

%% Implémentation des cas de charge équivalente dans le modèle

nodes = transpose(1:NNode);
fx = zeros(size(nodes));
fy = zeros(size(nodes));
mz = zeros(size(nodes));

for i = 1:NNode
    
    if P1_data(3*i-2) ~= 0
        fx(i) = P1_data(3*i-2);
    end
    if P1_data(3*i-1) ~=0
        fy(i) = P1_data(3*i-1);
    end
    if P1_data(3*i) ~= 0
        mz(i) = P1_data(3*i);
    end
end


new_load_case = struct('nodes', nodes, 'fx', fx, 'fy', fy, 'mz', mz);

Add_LoadCase('Test8.BeamZ', new_load_case,'Test8_copie.BeamZ','P1')



nodes = transpose(1:NNode);
fx = zeros(size(nodes));
fy = zeros(size(nodes));
mz = zeros(size(nodes));

P1_data = Psigma_data;

for i = 1:NNode
    
    if P1_data(3*i-2) ~= 0
        fx(i) = P1_data(3*i-2);
    end
    if P1_data(3*i-1) ~=0
        fy(i) = P1_data(3*i-1);
    end
    if P1_data(3*i) ~= 0
        mz(i) = P1_data(3*i);
    end
end

    
new_load_case = struct('nodes', nodes, 'fx', fx, 'fy', fy, 'mz', mz);

Add_LoadCase('Test8_copie.BeamZ', new_load_case,'Test8_copie.BeamZ','Psigma')