%=====================================================
%
% CS82 - 2012/2013
% 
% Calcul FEM de la déformée d'une colonne isostress
% par des éléments de type "barre" linéaire
%
%=====================================================


%===== Nettoyage de Matlab
close all;
clear all;
clc;


%===== Paramètres
Force=500E3;  %N
sigma=1E6;  %Pa

rho=2400;   %kg/m^3
E=25E9;     %Pa

g=10;       %m/s^2
L=10;       %m


%===== Trace de A(x)
nPts=1000;
x=zeros(nPts,1);
A=zeros(nPts,1);
R=zeros(nPts,1);%rayon du disque associé
for i=1:nPts
    xi=(i-1)*L/(nPts-1);
    x(i,1)=xi;
    A(i,1)=Aire(Force,sigma,rho,g,L,xi);
    R(i,1)=sqrt(A(i,1)/pi);
end

figure(1)
hold on
plot(x,R,'-b','linewidth',2);
grid;
xlabel('x (m)');
ylabel('R(x) (m)');
Title('Rayon de la colonne isotress de section circulaire'); 
hold off


%===== Creation de la matrice de rigidite globale
nNodes=40;
nElem=nNodes-1;
Le=L/nElem;

x=zeros(nNodes,1);
K_g=zeros(nNodes,nNodes);

for i=1:nNodes
    xe=(i-1)*Le;
    x(i,1)=xe;
end

for i=1:nElem
    xe=(i-1)*Le;
    Ae=Aire(Force,sigma,rho,g,L,xe);
    K_e=E*Ae/Le*[1,-1;-1,1];
    K_g(i,i)=       K_g(i,i)+       K_e(1,1);
    K_g(i,i+1)=     K_g(i,i+1)+     K_e(1,2);
    K_g(i+1,i)=     K_g(i+1,i)+     K_e(2,1);
    K_g(i+1,i+1)=   K_g(i+1,i+1)+   K_e(2,2);
end


%===== Creation du vecteur 2nd membre global
f_g=zeros(nNodes,1);

for i=1:nElem
    xe=(i-1)*Le;
    Ae=Aire(Force,sigma,rho,g,L,xe);
    f_e=0.5*rho*(-g)*Le*Ae*[1;1];
    f_g(i,1)=   f_g(i,1)+    f_e(1,1);
    f_g(i+1,1)= f_g(i+1,1)+  f_e(2,1);
end


%===== Imposition des conditions-limites
% CL cinématique u(x=0)=0
K_g(1,:)=zeros(1,nNodes);
K_g(:,1)=zeros(nNodes,1);
K_g(1,1)=1;
f_g(1,1)=0;

% CL en effort F(x=L)=-Force
f_g(nNodes,1)=f_g(nNodes,1)-Force;


%===== Résolution
u_g=zeros(nNodes,1);
u_g=K_g\f_g;

figure(2)
hold on
plot(x,u_g,'-go','linewidth',2);
xlabel('x (m)');
ylabel('u(x) (m)');
title('Deformee de la colonne');
grid;
hold off


%===== Post-traitement
stress=zeros(nElem,1);
for i=1:nElem
    epsilon=(u_g(i+1)-u_g(i))/Le;
    stress(i,1)=E*epsilon;
end

figure(3)
hold on
for i=1:nElem
    plot([x(i,1),x(i+1,1)],stress(i,1)*[1,1],'-ro','linewidth',2);
end
xlabel('x (m)');
ylabel('\sigma(x) (Pa)');
title('Contrainte normale de la colonne');
grid;
hold off