l = [0.8 1 1.2]; % Longueur  des éléments [m]
E = [210e9, 69e9, 80e9]; % Modules de Young [Pa]
A0 = 200e-6; % Section initiale [m^2]
lambda = 0.5; % Facteur exponentiel
x_mid = [0.4 1.3 2.4]; % milieu des éléments
A = A0 * exp(-lambda * x_mid); % Section au milieu de l'élément
F = [10000, -8000, 3000]; % Forces appliquées [N]

% Nœuds et éléments
nodes = [0 0.4 0.8 1.3 1.8 2.4 3.0];
num_nodes = length(nodes);
num_elements = 3;

% Préallocation des matrices
stiffness = zeros(num_nodes);
applied_forces = zeros(num_nodes,1);

%Fonctions de base quadratiqeus
syms xsi
he1 = symfun(xsi*(xsi-1)/2, xsi);
he2 = symfun(xsi*(xsi+1)/2, xsi);
he3 = symfun(1-xsi^2, xsi);

%Matrices élémentaires
He = [he1, he2, he3];
Be = diff(He,xsi);

stiffness_elem1 = (2*E(1)*A(1)/l(1))*int(transpose(Be)*Be,xsi,[-1 1]);
stiffness_elem2 = (2*E(2)*A(2)/l(2))*int(transpose(Be)*Be,xsi,[-1 1]);
stiffness_elem3 = (2*E(3)*A(3)/l(3))*int(transpose(Be)*Be,xsi,[-1 1]);

applied_forces_elem1 = [0,0,F(1)];
applied_forces_elem2 = [0,0,F(2)];
applied_forces_elem3 = [0,F(3),0];

%Assemblage 
function y = QuadraticBarAssembleMatrix(K,ke,i,j,m)
%QuadraticBarAssembleMatrix This function assembles the element stiffness
% matrix ke of the quadratic bar with nodes i, j
% and m into the global stiffness matrix K.
% This function returns the global stiffness
% matrix K after the element stiffness matrix
% k is assembled.
K(i,i) = K(i,i) + ke(1,1);
K(i,j) = K(i,j) + ke(1,2);
K(i,m) = K(i,m) + ke(1,3);
K(j,i) = K(j,i) + ke(2,1);
K(j,j) = K(j,j) + ke(2,2);
K(j,m) = K(j,m) + ke(2,3);
K(m,i) = K(m,i) + ke(3,1);
K(m,j) = K(m,j) + ke(3,2);
K(m,m) = K(m,m) + ke(3,3);
y = K;
end

function y = QuadraticBarAssembleVector(r,re,i,j,m)
%QuadraticBarAssembleVector This function assembles the element
% vector of applied forces re of the quadratic bar with nodes i, j
% and m into the global vector of applied forces r.
% This function returns the global vector of applied forces
% r after the element vector of applied forces
% re is assembled.
r(i) = r(i) + re(1);
r(j) = r(j) + re(2);
r(m) = r(m) + re(3);
y = r;
end

stiffness = QuadraticBarAssembleMatrix(stiffness,stiffness_elem1,1,3,2);
stiffness = QuadraticBarAssembleMatrix(stiffness,stiffness_elem2,3,5,4);
stiffness = QuadraticBarAssembleMatrix(stiffness,stiffness_elem3,5,7,6);

applied_forces= QuadraticBarAssembleVector(applied_forces,applied_forces_elem1,1,3,2);
applied_forces= QuadraticBarAssembleVector(applied_forces,applied_forces_elem2,3,5,4);
applied_forces= QuadraticBarAssembleVector(applied_forces,applied_forces_elem3,5,7,6);

% Conditions aux limites : encastrement à l'extrémité gauche (déplacement nul)
stiffness_free_dof = stiffness(2:end, 2:end);
force_free_dof = applied_forces(2:end);

% Résolution du système global
displacements_free_dof = stiffness_free_dof \ force_free_dof;
displacements = [0; displacements_free_dof];

% Affichage des résultats
syms x

for elem_index = 1:num_elements
    right_node_index = 2*elem_index+1;
    left_node_index = 2*elem_index-1;
    middle_node_index = 2*elem_index;
    
    right_node = nodes(right_node_index);
    left_node = nodes(left_node_index);
    
    xToxsiTransform = (2*x-(right_node+left_node))/(right_node-left_node);
    H = compose(He,xToxsiTransform);
    uh_local = H*[displacements(left_node_index); displacements(right_node_index); displacements(middle_node_index)];
    fplot(uh_local, [left_node right_node]);
    hold on
end
scatter(nodes, displacements, 'o');
xlabel('Position x [m]');
ylabel('Déplacement u [m]');
title('Déplacement axial');
hold off