
clear all
close all



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition des parametres utiles

alphaM = pi/2;
alphaVert  = pi/6;
alphaTd = pi/2 + alphaVert;
alphaTg = pi/2 - alphaVert;

T = -1000; % dans les tirants (< en traction)
Tmat = -2*T*cos(alphaVert); % dans le mat
Eacier = 2.1e11;
ro = 7800;
S1 = pi*0.005^2; % tirant

diamExterieurMat = 0.1;
diamInterieurMat = 0.1-2*0.005;

S2 = pi*((diamExterieurMat/2)^2-(diamInterieurMat/2)^2); % mat
J = pi*(diamExterieurMat^4-diamInterieurMat^4)/64;
L = 4 ;
Ltirant = L / cos(alphaVert);


%%%%%%%%%%%%%%%%%%%%%%%%
% Phase 1. Maillage
% (defini par l'enonce)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Phase 2. Formulation 
% (recuperation des resultats du cours)

syms Nsc Lc Ec Sc roSc   Nsp Lp Ep Sp roSp Ip   

% - construction des elements de reference de poutre Eb 2D
KpP=(-Nsp/(30*Lp))*[0 0 0 0 0 0; 0 36 3*Lp 0 -36 3*Lp ; 0 3*Lp 4*Lp^2 0 -3*Lp -Lp^2 ; 0 0 0 0 0 0 ; 0 -36 -3*Lp 0 36 -3*Lp ; 0 3*Lp -Lp^2 0 -3*Lp 4*Lp^2];
KeP=(Ep*Ip/Lp^3)*[Sp*Lp^2/Ip 0 0 -Sp*Lp^2/Ip 0 0 ; 0 12 6*Lp 0 -12 6*Lp ; 0 6*Lp 4*Lp^2 0 -6*Lp 2*Lp^2 ; -Sp*Lp^2/Ip 0 0 Sp*Lp^2/Ip 0  0;  0 -12 -6*Lp 0 12 -6*Lp ; 0 6*Lp 2*Lp^2 0 -6*Lp 4*Lp^2];
MP=(roSp*Lp/420)*[140 0 0 70 0 0 ; 0 156 22*Lp 0 54 -13*Lp ; 0 22*Lp 4*Lp^2 0 13*Lp -3*Lp^2 ; 70 0 0 140 0 0 ; 0 54 13*Lp 0 156 -22*Lp ; 0 -13*Lp -3*Lp^2 0 -22*Lp 4*Lp^2];

% - construction des elements de reference de cable tendu ( = barre 2D)
KpC=(-Nsc/Lc)*[0 0 0 0 ; 0 1 0 -1 ; 0 0 0 0 ; 0 -1 0 1];
KeC=(Ec*Sc/Lc)*[1 0 -1 0 ; 0 0 0 0 ; -1 0 1 0 ; 0 0 0 0];
MC=(roSc*Lc/6)*[2 0 1 0 ; 0 2 0 1 ; 1 0 2 0 ; 0 1 0 2];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Phase 3. Assemblage des sous-domaines

% - reformulation des matrices elementaires dans le repere global

% ----> rotation du mat
%RMat = [cos(alphaM) sin(alphaM) 0 ; -sin(alphaM) cos(alphaM) 0 ; 0 0 1];
RMat = [0 1 0 ; -1 0 0 ; 0 0 1];
TMat = [RMat zeros(3,3) ; zeros(3,3) RMat];
KpMat = TMat' * KpP *TMat; 
KeMat = TMat' * KeP *TMat; 
MMat = TMat' * MP *TMat; 

% ----> affectation des parametres du mat 
KpMat = subs(KpMat,{Nsp, Lp, Ep, Ip, Sp, roSp}, {Tmat, L, Eacier, J, S2, ro*S2});
KeMat = subs(KeMat,{Nsp, Lp, Ep, Ip, Sp, roSp}, {Tmat, L, Eacier, J, S2, ro*S2});
MMat = subs(MMat,{Nsp, Lp, Ep, Ip, Sp, roSp}, {Tmat, L, Eacier, J, S2, ro*S2});


% ----> rotation du tirant droit
RTd = [cos(alphaTd) sin(alphaTd) ; -sin(alphaTd) cos(alphaTd)];
TTd = [RTd zeros(2,2) ; zeros(2,2) RTd];
KpTd = TTd' * KpC *TTd; 
KeTd = TTd' * KeC *TTd; 
MTd = TTd' * MC *TTd; 

% ----> rotation du tirant gauche
RTg = [cos(alphaTg) sin(alphaTg) ; -sin(alphaTg) cos(alphaTg)];
TTg = [RTg zeros(2,2) ; zeros(2,2) RTg];
KpTg = TTg' * KpC *TTg; 
KeTg = TTg' * KeC *TTg; 
MTg = TTg' * MC *TTg; 

% ----> affectation des parametres des tirants
KpTd = subs( KpTd, {Nsc, Lc, Ec, Sc, roSc},{T, Ltirant, Eacier, S1, ro*S1});
KeTd = subs( KpTd, {Nsc, Lc, Ec, Sc, roSc},{T, Ltirant, Eacier, S1, ro*S1});
MTd = subs( KpTd, {Nsc, Lc, Ec, Sc, roSc},{T, Ltirant, Eacier, S1, ro*S1});
KpTg = subs( KpTd, {Nsc, Lc, Ec, Sc, roSc},{T, Ltirant, Eacier, S1, ro*S1});
KeTg = subs( KpTd, {Nsc, Lc, Ec, Sc, roSc},{T, Ltirant, Eacier, S1, ro*S1});
MTg = subs( KpTd, {Nsc, Lc, Ec, Sc, roSc},{T, Ltirant, Eacier, S1, ro*S1});


% - expression de la compatibilite des deplacements nodaux (matrice de localisation, connectivite)

% ----> placement des ddl du mat parmi ceux de la structure
Bmat = [eye(6) zeros(6, 4)];

% ----> placement des ddl du tirant droit parmi ceux de la structure
Btd = [0 0 0 0 0 0 1 0 0 0 ; 0 0 0 0 0 0 0 1 0 0 ; 0 0 0 1 0 0 0 0 0 0 ; 0 0 0 0 1 0 0 0 0 0];

% ----> placement des ddl du tirant gauche parmi ceux de la structure
Btg = [0 0 0 0 0 0 0 0 1 0 ; 0 0 0 0 0 0 0 0 0 1 ; 0 0 0 1 0 0 0 0 0 0 ; 0 0 0 0 1 0 0 0 0 0];

% ----> matrices assemblees, globales, de la structure
Kpglobale = Bmat'*KpMat*Bmat + Btd'*KpTd*Btd + Btg'*KpTg*Btg; 
Keglobale = Bmat'*KeMat*Bmat + Btd'*KeTd*Btd + Btg'*KeTg*Btg; 
Mglobale = Bmat'*MMat*Bmat + Btd'*MTd*Btd + Btg'*MTg*Btg; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Phase 4. Prise en compte des conditions aux limites cinematiques

% - prise en compte des CL cinematiques
ddlAexclure = [1 2 3 5 7 8 9 10] ; % liste des ddl a exclure
Kpglobale(ddlAexclure,:)=[]; Kpglobale(:,ddlAexclure)=[] 
Keglobale(ddlAexclure,:)=[]; Keglobale(:,ddlAexclure)=[]
Mglobale(ddlAexclure,:)=[]; Mglobale(:,ddlAexclure)=[]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Phase 5. Resolution (calcul des modes)

[V , D] = eig(Keglobale+Kpglobale,Mglobale);
pulsations = sqrt(diag(D))
V(:,1);
V(:,2);

xmat = linspace(0,1,20);
N5 = 3*xmat.^2-2*xmat.^3;
N6 = -L*(xmat.^2 -xmat.^3);

mode1 = V(1,1)*N5 +V(2,1)*N6;
mode2 = V(1,2)*N5 +V(2,2)*N6;

figure
plot(mode1)
figure
plot(mode2)





