%Program to analyze trusses with supports and forces specified.
%Specified support displacements are not allowed.
clear

%*** INPUT STARTS HERE ***
%Specify structure properties
nnodes=4; %specify the number of nodes
nmem=5; %specify the number of members

%Define nodal coordinates in order starting with node one
x(1)=4; y(1)=3;
x(2)=4; y(2)=0;
x(3)=8; y(3)=0;
x(4)=0; y(4)=0;

%Define member connectivity start with member one
% Start node        end node
mconn(1,1)=4; mconn(1,2)=1;
mconn(2,1)=1; mconn(2,2)=3;
mconn(3,1)=2; mconn(3,2)=3;
mconn(4,1)=4; mconn(4,2)=2;
mconn(5,1)=2; mconn(5,2)=1;

%input properties for each member

%For now E is constant.
E=200*10^6;

%A for each member
A(1)=0.0015;
A(2)=0.0015;
A(3)=0.0015;
A(4)=0.0015;
A(5)=0.0015;

%Specify supports
nsupp=3; %number of degrees of freedom with supports
sup(1)=7;%for the 1 to nsupp indicate which dof are supports
sup(2)=8;
sup(3)=6;
%sup(4)=6;
%sup(5)=7;
%sup(6)=8;

%Specify forces
%zero all the forces
for i=1:nnodes*2;
    F(i)=0.0;
end;    
nforce=5; %number of degrees of freedom with forces given
F(1)=0; %specify force at each degree of freedom where the external force is known, F(dof)
F(2)=0;
F(3)=0;
F(4)=-10;
F(5)=0

%***INPUT ENDS HERE ***

%Calculate the truss member lengths, and angles
for i=1:nmem;
    dx=x(mconn(i,2))-x(mconn(i,1));
    dy=y(mconn(i,2))-y(mconn(i,1));
    L(i)=sqrt(dx^2+dy^2);
    c(i)=dx/L(i);
    s(i)=dy/L(i);
end;

%For each element create the global stiffness matrix and put it into the
%global stiffness matrix.

%zero the global stiffness matrix
for i=1:nnodes*2;
    for j=1:nnodes*2;
        kg(i,j)=0.0;
    end;
end;    

%Create each element global stiffness matrix as transpose[T][k][T]
for m=1:nmem; %loop over each element
    
    %zero k and T
    for i=1:4;
        for j=1:4;
            k(i,j)=0.0;
            T(i,j)=0.0;
        end;    
    end;
    
    %create T
    T(1,1)=c(m) ; T(1,2)=s(m); T(1,3)=0;     T(1,4)=0;
    T(2,1)=-s(m); T(2,2)=c(m); T(2,3)=0;     T(2,4)=0;
    T(3,1)=0;     T(3,2)=0;    T(3,3)=c(m);  T(3,4)=s(m);
    T(4,1)=0;     T(4,2)=0;    T(4,3)=-s(m); T(4,4)=c(m);
    
    %create k    
    k(1,1)=1.0; k(1,3)=-1.0; k(3,1)=-1.0; k(3,3)=1.0;
    k=A(m)*E/L(m)*k;
    
    %transform k into the element global stiffness matrix
    k=T'*k*T;
    
    %put each element k into the global stiffness matrix kg
    for i=1:2;
        for j=1:2;
            kg(mconn(m,i)*2-1,mconn(m,j)*2-1)=kg(mconn(m,i)*2-1,mconn(m,j)*2-1)+k(i*2-1,j*2-1);
            kg(mconn(m,i)*2,mconn(m,j)*2)=kg(mconn(m,i)*2,mconn(m,j)*2)+k(i*2,j*2);
            kg(mconn(m,i)*2-1,mconn(m,j)*2)=kg(mconn(m,i)*2-1,mconn(m,j)*2)+k(i*2-1,j*2);
            kg(mconn(m,i)*2,mconn(m,j)*2-1)=kg(mconn(m,i)*2,mconn(m,j)*2-1)+k(i*2,j*2-1);
        end;    
    end; %now kg should be complete
    
end; %end loop over all the elements    

%Put kg into kgs and modify kgs to solve for the unknown displacements
%we save kg so that we can calculate the reactions also.
kgs=kg;

%modify kgs based on the support conditions.  Zero the entire row for each
%dof sup(i) except for entry kgs(sup(i),sup(i)) set it equal to 1.  When we solve
%for the displacements we have essentially set them equal to zero since
%they are supports.
for i=1:nsupp; %loop through all the specified supports
    for j=1:nnodes*2;
            kgs(sup(i),j)=0.0; %zero the whole row
    end;
    kgs(sup(i),sup(i))=1.0; %put a 1 in position kgs(sup(i),sup(i))
end;
kgs
%solve for the global nodal displacements
display('The global nodal displacements are:')
U=inv(kgs)*F'

%Solve for the global external forces, here we use kg that we set aside
%above.
display('The global external forces are:')
F=kg*U

%Calculate the forces for each member
for m=1:nmem; %loop over each element
    
    %zero k and T
    for i=1:4;
        for j=1:4;
            k(i,j)=0.0;
            T(i,j)=0.0;
        end;    
    end;
    
    %create T
    T(1,1)=c(m) ; T(1,2)=s(m); T(1,3)=0;     T(1,4)=0;
    T(2,1)=-s(m); T(2,2)=c(m); T(2,3)=0;     T(2,4)=0;
    T(3,1)=0;     T(3,2)=0;    T(3,3)=c(m);  T(3,4)=s(m);
    T(4,1)=0;     T(4,2)=0;    T(4,3)=-s(m); T(4,4)=c(m);
    
    %create k    
    k(1,1)=1.0; k(1,3)=-1.0; k(3,1)=-1.0; k(3,3)=1.0;
    k=A(m)*E/L(m)*k;
    
    %Extract the global displacements associated with element m
    u(1)=U(mconn(m,1)*2-1);
    u(2)=U(mconn(m,1)*2);
    u(3)=U(mconn(m,2)*2-1); 
    u(4)=U(mconn(m,2)*2);
    
    %calculate the local element end forces
    f=k*T*u';
    
    %save the important end forces
    fm(m)=f(3); %positive values mean truss member in tension, negative means compression
end;%end loop over each element

display('The member forces are:')
fm'

%plot the truss on a graph
plot(0,0,'w')
hold
%loop over all the elements and plot each element, one at a time
magnify=2000; %magnification factor for deflections
for m=1:nmem;
    xm(1)=x(mconn(m,1));
    ym(1)=y(mconn(m,1));
    xm(2)=x(mconn(m,2));
    ym(2)=y(mconn(m,2));
    xu(1)=U(mconn(m,1)*2-1)*magnify;
    yu(1)=U(mconn(m,1)*2)*magnify;
    xu(2)=U(mconn(m,2)*2-1)*magnify;
    yu(2)=U(mconn(m,2)*2)*magnify;
    if fm(m) < 0.0;
    plot(xm,ym,'k'); %undisplaced truss members in compression will be blue
    else;
    plot(xm,ym,'b'); %undisplaced truss members in tension will be red 
    end;
    plot(xm+xu,ym+yu,'r--');  %displaced truss members
end;
xlabel('xcoordinates')
ylabel('ycoordinates')
title('Undeflected(solid) and Magnified Deflected(dashed) Truss')
hold


