%=====================================================
%
% CS82 - 2012/2013
% 
% Calcul FEM de la déformée d'un treillis constitué
% d'éléments de type "barre" 2D linéaire
%
%=====================================================


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


%===== Paramètres
F=-10000;   %N
E=210E9;%Pa
A=1E-4; %m^2
L=1;    %m


%===== Maillage
% %Maillage exo analytique
% coor=[0,L;L,L;0,0];
% conn=[1,2;3,2];
% nodesBloques=[1;3];
% iforce=4;%indice dans le vecteur global

%Treillis 10 barres
coor=[0,0;L,0;2*L,0;0,L;L,L;2*L,L];
conn=[1,2;2,3;1,5;2,4;2,6;3,5;4,5;5,6;3,6;2,5];
nodesBloques=[1;4];
iforce=6;

[nNodes,temp]=size(coor);
[nElem,temp]=size(conn);

figure(1)
hold on
for i=1:nNodes
    plot(coor(i,1),coor(i,2),'bo','linewidth',2);
    text(coor(i,1),coor(i,2)+5E-2,num2str(i),'color',[0,0,0]);
end
for i=1:nElem
    plot(coor(conn(i,1:2),1),coor(conn(i,1:2),2),'-b','linewidth',2);
    text(sum(coor(conn(i,1:2),1))/2,sum(coor(conn(i,1:2),2))/2-5E-2,num2str(i),'color',[0,0,1]);
end
grid;
axis equal;
xlabel('x (m)');
ylabel('y (m)');
Title('Maillage du treillis'); 
hold off


%===== Creation des donnees intermediaires
tabL=zeros(nElem,1);
tabAngle=zeros(nElem,1);
for i=1:nElem
    x1=coor(conn(i,1),1);
    y1=coor(conn(i,1),2);
    x2=coor(conn(i,2),1);
    y2=coor(conn(i,2),2);
    
    tabL(i,1)=sqrt( (x2-x1)*(x2-x1)+(y2-y1)*(y2-y1) );
    tabAngle(i,1)=atan2(y2-y1,x2-x1);
end


%===== Creation de la matrice de rigidite globale
K=zeros(2*nNodes,2*nNodes);

for i=1:nElem
    Ke=E*A/tabL(i,1)*[1,-1;-1,1];
    c=cos(tabAngle(i,1));
    s=sin(tabAngle(i,1));
    Te=[c,s,0,0;0,0,c,s];
    Ke=Te'*Ke*Te;
    i1=conn(i,1);
    i2=conn(i,2);
    indices=[2*(i1-1)+1;2*(i1-1)+2;2*(i2-1)+1;2*(i2-1)+2];
    for j=1:4
        ij=indices(j,1);
        for k=1:4
            ik=indices(k,1);
            K(ij,ik)=K(ij,ik)+Ke(j,k);
        end
    end
end


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


%===== Imposition des conditions-limites
% CL cinématique
for i=1:length(nodesBloques)
    for j=1:2
        ii=2*(nodesBloques(i,1)-1)+j;
        K(ii,:)=zeros(1,2*nNodes);
        K(:,ii)=zeros(2*nNodes,1);
        K(ii,ii)=1;
    end
end

% CL en effort
f(iforce,1)=F;%au cas par cas


%===== Résolution
u=zeros(nNodes,1);
u=K\f;


%===== Affichage de la deformee
amplification=100;

figure(2)
hold on
for i=1:nNodes
    plot(coor(i,1),coor(i,2),'bo','linewidth',2);
    text(coor(i,1),coor(i,2)+5E-2,num2str(i),'color',[0,0,0]);
end
for i=1:nElem
    plot(coor(conn(i,1:2),1),coor(conn(i,1:2),2),'-b','linewidth',2);
    text(sum(coor(conn(i,1:2),1))/2,sum(coor(conn(i,1:2),2))/2-5E-2,num2str(i),'color',[0,0,1]);
end

for i=1:nNodes
    coor(i,1)=coor(i,1)+amplification*u(2*(i-1)+1,1);
    coor(i,2)=coor(i,2)+amplification*u(2*(i-1)+2,1);
end

for i=1:nNodes
    plot(coor(i,1),coor(i,2),'ro','linewidth',2);
    text(coor(i,1),coor(i,2)+5E-2,num2str(i),'color',[1,0,0]);
end
for i=1:nElem
    plot(coor(conn(i,1:2),1),coor(conn(i,1:2),2),'-r','linewidth',2);
    text(sum(coor(conn(i,1:2),1))/2,sum(coor(conn(i,1:2),2))/2-5E-2,num2str(i),'color',[1,0,0]);
end

grid;
axis equal;
xlabel('x (m)');
ylabel('y (m)');
Title('Deformee du treillis'); 
hold off