%% 
% Fonction qui sert à ressort les résultats grâce à la méthode des éléments
% finis.

function [K0, Ksigma, phi0, L_FEM, Re, Model,ddlLibres] = calcul_FEM(Model,listDDLbloque)


NNode = length(Model.Nodes.XNOD);
NElem = length(Model.Elements.ELEMNOA);

% K0 - matrice de raideur assemblée pour modèle non perturbé
[K0,Ksigma,~,~,~,~,~, Ke,~,Re,Model.Elements.ELEMLEN,Model.Elements.ELEMDOF] = ...
    Assemble(NNode,NElem,Model.Nodes,Model.Elements, Model.Geometries, Model.Materials, Model.Springs);

[P,Pkip] = AssembleLoads(NNode,NElem,Model.Nodes,Model.Elements,Re,Model.Loads,Model.Geometries,Model.Materials,Model.Springs);
P = P(:,2); % on conserve uniquement le 2eme cas de charge
P_all = P;
Pkip = Pkip(:,2);

ddlLibres = setdiff(1:3*NNode,listDDLbloque);
K00 = K0;
K0 = K0(ddlLibres,ddlLibres);
P = P(ddlLibres,1);
x = K0\P;

X = zeros(3*NNode,1);
X(ddlLibres)=x;

for iel=1:NElem
    xloc = X(Model.Elements.ELEMDOF(:,iel),1);
    Fint = Re(:,:,iel)*Ke(:,:,iel)*xloc;
    Model.Elements.N_STAB(iel) = -Fint(1); % save axial force for stability analysis
end

%Recalcul en prenant en compte les forces

[K0,Ksigma] = ...
    Assemble(NNode,NElem,Model.Nodes,Model.Elements, Model.Geometries, Model.Materials, Model.Springs);

K0 = K0(ddlLibres,ddlLibres);
Ksigma1 = Ksigma;
Ksigma = Ksigma(ddlLibres,ddlLibres);

[V0,Lambda0] = eig(K0,Ksigma);

phi0 = zeros(3*NNode,size(V0,2));
phi0(ddlLibres,:) = V0;

[Lambda0,ordre] = sort(diag(Lambda0));
phi0 = phi0(:,ordre);
ii = find(Lambda0>0);
Lambda0 = Lambda0(ii);
L_FEM = Lambda0(1);
phi0 = phi0(:,ii);
phi0 = phi0(:,1);
phi0 = phi0/max(abs([phi0(1:3:end); phi0(2:3:end)]));

K0 = K00;
Ksigma = Ksigma1;

end