%Program to analyze 2D - frames 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=3; %specify the number of members

%Define nodal coordinates in order starting with node one, node i must have
%degrees of freedom 3*i-2, 3*i-1, & 3*i associated with it.
x(1)=0; y(1)=144;
x(2)=96; y(2)=144;
x(3)=0; y(3)=0;
x(4)=96; y(4)=0;

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

%input properties for each member

%For now E is constant.
E=29*10^3;

%A and I for each member
A(1)=8; I(1)=450; %member 1
A(2)=8; I(2)=450; %member 2
A(3)=8; I(3)=450; %...

%Specify supports
nsupp=6; %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)=9;
sup(4)=10;
sup(5)=11;
sup(6)=12;

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

%***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*3;
    for j=1:nnodes*3;
        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:6;
        for j=1:6;
            k(i,j)=0.0;
            T(i,j)=0.0;
        end;    
    end;
    
    %create T, the transformation matrix
    T(1,1)=c(m) ; T(1,2)=s(m); T(1,3)=0; T(1,4)=0;     T(1,5)=0;    T(1,6)=0;
    T(2,1)=-s(m); T(2,2)=c(m); T(2,3)=0; T(2,4)=0;     T(2,5)=0;    T(2,6)=0;
    T(3,1)=0;     T(3,2)=0;    T(3,3)=1; T(3,4)=0;     T(3,5)=0;    T(3,6)=0;
    T(4,1)=0;     T(4,2)=0;    T(4,3)=0; T(4,4)=c(m);  T(4,5)=s(m); T(4,6)=0;
    T(5,1)=0;     T(5,2)=0;    T(5,3)=0; T(5,4)=-s(m); T(5,5)=c(m); T(5,6)=0;
    T(6,1)=0;     T(6,2)=0;    T(6,3)=0; T(6,4)=0;     T(6,5)=0;    T(6,6)=1;
    
    %create k, the local element stiffness matrix    
    k(1,1)=A(m)*E/L(m); k(1,2)=0; k(1,3)=0; k(1,4)=-A(m)*E/L(m); k(1,5)=0; k(1,6)=0;
    k(2,1)=0; k(2,2)=12*E*I(m)/L(m)^3; k(2,3)=6*E*I(m)/L(m)^2; k(2,4)=0; k(2,5)=-12*E*I(m)/L(m)^3; 
    k(2,6)=k(2,3);
    k(3,1)=0; k(3,2)=k(2,3); k(3,3)=4*E*I(m)/L(m); k(3,4)=0; k(3,5)=-k(2,3); k(3,6)=2*E*I(m)/L(m);
    k(4,1)=-A(m)*E/L(m); k(4,2)=0; k(4,3)=0; k(4,4)=A(m)*E/L(m); k(4,5)=0; k(4,6)=0;    
    k(5,1)=0; k(5,2)=k(2,5); k(5,3)=k(3,5); k(5,4)=0; k(5,5)=k(2,2); k(5,6)=-k(2,3);
    k(6,1)=0; k(6,2)=k(2,6); k(6,3)=k(3,6); k(6,4)=k(4,6); k(6,5)=k(5,6); k(6,6)=k(3,3);
    
    %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)*3-2,mconn(m,j)*3-2)=kg(mconn(m,i)*3-2,mconn(m,j)*3-2)+k(i*3-2,j*3-2);
            kg(mconn(m,i)*3-1,mconn(m,j)*3-1)=kg(mconn(m,i)*3-1,mconn(m,j)*3-1)+k(i*3-1,j*3-1);
            kg(mconn(m,i)*3,mconn(m,j)*3)=kg(mconn(m,i)*3,mconn(m,j)*3)+k(i*3,j*3);
            
            kg(mconn(m,i)*3-2,mconn(m,j)*3-1)=kg(mconn(m,i)*3-2,mconn(m,j)*3-1)+k(i*3-2,j*3-1);
            kg(mconn(m,i)*3-2,mconn(m,j)*3)=kg(mconn(m,i)*3-2,mconn(m,j)*3)+k(i*3-2,j*3);
            
            kg(mconn(m,i)*3-1,mconn(m,j)*3-2)=kg(mconn(m,i)*3-1,mconn(m,j)*3-2)+k(i*3-1,j*3-2);
            kg(mconn(m,i)*3-1,mconn(m,j)*3)=kg(mconn(m,i)*3-1,mconn(m,j)*3)+k(i*3-1,j*3);
            
            kg(mconn(m,i)*3,mconn(m,j)*3-2)=kg(mconn(m,i)*3,mconn(m,j)*3-2)+k(i*3,j*3-2);        
            kg(mconn(m,i)*3,mconn(m,j)*3-1)=kg(mconn(m,i)*3,mconn(m,j)*3-1)+k(i*3,j*3-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*3;
            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;
display('The global stiffness matrix')
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:6;
        for j=1:6;
            k(i,j)=0.0;
            T(i,j)=0.0;
        end;    
    end;
    
    %create T, the transformation matrix
    T(1,1)=c(m) ; T(1,2)=s(m); T(1,3)=0; T(1,4)=0;     T(1,5)=0;    T(1,6)=0;
    T(2,1)=-s(m); T(2,2)=c(m); T(2,3)=0; T(2,4)=0;     T(2,5)=0;    T(2,6)=0;
    T(3,1)=0;     T(3,2)=0;    T(3,3)=1; T(3,4)=0;     T(3,5)=0;    T(3,6)=0;
    T(4,1)=0;     T(4,2)=0;    T(4,3)=0; T(4,4)=c(m);  T(4,5)=s(m); T(4,6)=0;
    T(5,1)=0;     T(5,2)=0;    T(5,3)=0; T(5,4)=-s(m); T(5,5)=c(m); T(5,6)=0;
    T(6,1)=0;     T(6,2)=0;    T(6,3)=0; T(6,4)=0;     T(6,5)=0;    T(6,6)=1;
    
    %create k, the element stiffness matrix in local coordinates    
    k(1,1)=A(m)*E/L(m); k(1,2)=0; k(1,3)=0; k(1,4)=-A(m)*E/L(m); k(1,5)=0; k(1,6)=0;
    k(2,1)=0; k(2,2)=12*E*I(m)/L(m)^3; k(2,3)=6*E*I(m)/L(m)^2; k(2,4)=0; k(2,5)=-12*E*I(m)/L(m)^3; 
    k(2,6)=k(2,3);
    k(3,1)=0; k(3,2)=k(2,3); k(3,3)=4*E*I(m)/L(m); k(3,4)=0; k(3,5)=-k(2,3); k(3,6)=2*E*I(m)/L(m);
    k(4,1)=-A(m)*E/L(m); k(4,2)=0; k(4,3)=0; k(4,4)=A(m)*E/L(m); k(4,5)=0; k(4,6)=0;    
    k(5,1)=0; k(5,2)=k(2,5); k(5,3)=k(3,5); k(5,4)=0; k(5,5)=k(2,2); k(5,6)=-k(2,3);
    k(6,1)=0; k(6,2)=k(2,6); k(6,3)=k(3,6); k(6,4)=k(4,6); k(6,5)=k(5,6); k(6,6)=k(3,3);
    
    %Extract the global displacements associated with element m
    u(1)=U(mconn(m,1)*3-2);
    u(2)=U(mconn(m,1)*3-1);
    u(3)=U(mconn(m,1)*3); 
    u(4)=U(mconn(m,2)*3-2);
    u(5)=U(mconn(m,2)*3-1);
    u(6)=U(mconn(m,2)*3);
    
    %calculate the local element end forces
    display ('Local End Forces for element number:')
    m
    f=k*T*u'
        
end;%end loop over each element

%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=200; %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)*3-2)*magnify;
    yu(1)=U(mconn(m,1)*3-1)*magnify;
    xu(2)=U(mconn(m,2)*3-2)*magnify;
    yu(2)=U(mconn(m,2)*3-1)*magnify;
    plot(xm,ym,'b'); %undisplaced frame members solid blue
    %Note that frame member displacements plotted here are not the true
    %displacements because the flexural displacements are not calculated.
    %Here we only take into account translation displacements at each node.
    plot(xm+xu,ym+yu,'r--');  %displaced frame members dashed red
end;
xlabel('xcoordinates')
ylabel('ycoordinates')
title('Undeflected(solid) and Magnified Deflected(dashed) Frame')
hold


