% File Name: Truss_Main
% Code to generate stiffness matrices & calculate displacements,strains,
% stresses etc for the truss structure given
clc;
clear all

%--------------Input the known------------------------%

a=input('Area of the truss element (in m^2):');                     
E=input('Young-s Modulus of the truss element (in Pa):');            
l=input('Length of the shortest truss element (in m):');                  

F6 = input('Force on node 3 in vertical direction (in N):');  
F9 = input('Force on node 5 in horizontal direction (in N):'); 
F11= input('Force on node 6 in horizontal direction (in N):');
F12= input('Force on node 3 in vertical direction (in N):'); 

%--------------Calculations------------------------%
k1=a*E/l;                  %stiffness constant = AE/L of the short trusses which are either horizontal or vertical links
k2=a*E/(sqrt(2)*l);         %stiffness constant = AE/L of the long trusses which are diagonal

%building different element stiffness matrices by calling the routine
%"stiff" - the parameters sent are angle in degrees, stiffness constant and
%the indices of the displacements for each element. The function gives a
%12x12 matrix with each omponent of the elemental stiffness matrix places
%at appropriate location.

K1g=stiff(0,k1,1,2,5,6);     %element 1, short truss, horizontal, displacements U1,U2,U5 and U6 - refer to notes for the picture/node numbering
K2g=stiff(45,k2,1,2,7,8);    %element 2, long truss, at 45 degree, displacements U1,U2,U7 and U8
K3g=stiff(0,k1,3,4,7,8);     %element 3, short truss, horizontal, displacements U3,U4,U7 and U8
K4g=stiff(90,k1,5,6,7,8);    %element 4, short truss, vertical, displacements U5,U6,U7 and U8
K5g=stiff(0,k1,5,6,9,10);    %element 5, short truss, horizontal, displacements U5,U6,U9 and U10
K6g=stiff(135,k2,9,10,7,8);  %element 6, long truss, at 135 degree, displacements U9,U10,U7 and U8
K7g=stiff(0,k1,7,8,11,12);   %element 7, short truss, horizontal, displacements U7,U8,U11 and U12
K8g=stiff(90,k1,9,10,11,12); %element 8, short truss, vertical, displacements U9,U10,U11 and U12

% assembly of the global stiffness matrix by adding different elemenal
% "global" stiffness matrices
KG=K1g+K2g+K3g+K4g+K5g+K6g+K7g+K8g;

% -------------Solution for unknown displacements and forces------------%
KG_act=KG([5:12],[5:12]);           % active component of the global stiffness matrix
F_act=[0;F6;0;0;F9;0;F11;F12];   % Forces at the nodes - refer to the notes for the picture/node numbering
U_act=inv(KG_act)*F_act;              % calculation of unknown displacements 
U_tot=[0;0;0;0;U_act];               % entire displacement vector; placing of the known displacements 
F_tot=KG*U_tot;                      % calculation of all forces - including constraints and given - this is a bit less cumbersome to write

%For display as the results
ndisp=[U_tot(1),U_tot(2);U_tot(3),U_tot(4);U_tot(5),U_tot(6);U_tot(7),U_tot(8);U_tot(9),U_tot(10);U_tot(11),U_tot(12)];  %Ux and Uy for one node in one row
nforce=[F_tot(1),F_tot(2);F_tot(3),F_tot(4);F_tot(5),F_tot(6);F_tot(7),F_tot(8);F_tot(9),F_tot(10);F_tot(11),F_tot(12)]; %Fx and Fy for one node in one row

disp('   Node no    Ux(mm)    Uy(mm)    Fx(kN)    Fy(kN)');
noderesult=[(1:6)',ndisp*10^3,nforce/(10^3)];
disp(noderesult);

% -------------Calculation of stress and strain------------%
%Function StressStrain is invoked to achieve stress and strain for each
%element. The function takes the link number, angle wrt horizontal in
%degrees, diplacement at each nodes, youngs modulus and length of the link
[eps1, sig1]=StressStrain(0,U_tot([1,2,5,6],:),E,l);            %element 1, short truss, horizontal, displacements U1,U2,U5 and U6 
[eps2, sig2]=StressStrain(45,U_tot([1,2,7,8],:),E,sqrt(2)*l);   %element 2, long truss, at 45 degree, displacements U1,U2,U7 and U8 
[eps3, sig3]=StressStrain(0,U_tot([3,4,7,8],:),E,l);            %element 3, short truss, horizontal, displacements U3,U4,U7 and U8
[eps4, sig4]=StressStrain(90,U_tot([5,6,7,8],:),E,l);           %element 4, short truss, vertical, displacements U5,U6,U7 and U8
[eps5, sig5]=StressStrain(0,U_tot([5,6,9,10],:),E,l);           %element 5, short truss, horizontal, displacements U5,U6,U9 and U10
[eps6, sig6]=StressStrain(135,U_tot([9,10,7,8],:),E,sqrt(2)*l); %element 6, long truss, at 135 degree, displacements U9,U10,U7 and U8
[eps7, sig7]=StressStrain(0,U_tot([7,8,11,12],:),E,l);          %element 7, short truss, horizontal, displacements U7,U8,U11 and U12
[eps8, sig8]=StressStrain(90,U_tot([9,10,11,12],:),E,l);        %element 8, short truss, vertical, displacements U9,U10,U11 and U12

elm_strain=[eps1; eps2; eps3; eps4; eps5; eps6; eps7; eps8];    %strain column, strain of each link is a row
elm_stress=[sig1; sig2; sig3; sig4; sig5; sig6; sig7; sig8];    %stress column, stress of each link is a row
linkresults=[single(1:8)',elm_strain, elm_stress/10^6];

disp('Link no       Strain  Stress(MPa)');
disp(linkresults);