clear all

%%
% Ce code est le code général pour calculer n'importe quel modèle sorti de
% BeamZ. Nous pouvons comparer le résultat du multiplicateur critique
% approché par notre méthode au résultat exact calculé par la méthode des
% éléments finis pour différentes raideurs de ressort. Il faut néanmoins
% s'assurer que nous initialisons bien les DDL à bloquer pour que ça colle
% avec le modèle .BeamZ.


% Charger les données d'un fichier .mat
filename = 'Test4_symetric.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];

% 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];

% Initialiser les vecteurs 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);
e = zeros(length(K_ressort),1);

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

% Calcul de certaines constantes pour la suite
A = Model.Geometries.A;
E = Model.Materials(1).Young;
EA = E*A;
I = E;

% 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);
phi00 = phi0;
L0 = L_FEM;

% Boucle pour calculer le lambda approximatif en fonction de K
for i = 1:length(K_ressort)
    
    % Modifier le modèle avec la nouvelle constante de raideur pour les ressorts
    Model.Springs(3).Value = K_ressort(i);
    EI = Model.Geometries.I;
    
    
    % Calcul de Lambda1
    DH0 = 0;
    DHsigma = 0;
    Hsigma = 0;
    phi0 = phi00;

    % Parcourir tous les éléments pour calculer Lambda1
    for j = 1:NElem
        %Données
        phi0 = phi00;
        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
            alpha1 = 1/(K_ressort(i)*Le);
        end
        if Model.Elements.ELEMSP2(j) == 3
            alpha2 = 1/(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 = DH0 - 4*EI/Le*(alpha1*(2*phi1 + phi2 - 3*psi)^2 + alpha2*(phi1 + 2*phi2 - 3*psi)^2);

            %Calcul de DHsigma
            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 = Hsigma + P*Le/15 * (2*phi1^2 - phi1*phi2 + 2*phi2^2 - 3*(phi1+phi2)*psi + 18*psi^2);

    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;



end





%% 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_L_FEM, 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 8);
grid on;
legend('$\lambda_{approx}$', '$\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([2 5]);

% Ajout des valeurs de l'erreur en pourcentage pour chaque point
for i = 1:length(K_ressort)
    text(K_ressort(i), max(vect_sum(i), vect_L_FEM(i))+0.2*(min(vect_sum(1), vect_L_FEM(1))-max(vect_sum(1), vect_L_FEM(1))), sprintf('%.2f %%', err_pct(i)), 'HorizontalAlignment', 'center', 'FontSize', 9, 'Color', 'black');
end

% Ajout de la ligne et de la valeur de l'erreur maximale
text(K_ressort(idx)+1, min(vect_sum(idx), vect_L_FEM(idx))+0.1*(max(vect_sum(idx), vect_L_FEM(idx))-min(vect_sum(idx), vect_L_FEM(idx))), sprintf('Erreur max : %.2f %%', err_pct(idx)));