%Program to analyze 2D - frames with supports and forces specified.
%Specified support displacements are not allowed.
%Supports distributed loads and point loads between the nodes.
clear

%*** INPUT STARTS HERE ***
%Specify structure properties
nnodes=3; %specify the number of nodes
nmem=2; %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)=96;
x(2)=144; y(2)=96;
x(3)=144; y(3)=0;

%Define member connectivity start with member one
% Start node        end node
mconn(1,1)=1; mconn(1,2)=2; %member 1
mconn(2,1)=2; mconn(2,2)=3; %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)=10; I(1)=600; %member 1
A(2)=10; I(2)=600; %member 2
%A(3)=8; I(3)=450; %...

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

%(no input here)Calculate the frame 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;

%Specify joint forces
%zero all the joint forces
for i=1:nnodes*3;
    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)=0;
F(5)=0;
%F(6)=0;

%Specify forces on members in local coordinates
%A point load(at distance da from node 1) between member nodes and a uniformly distributed load
%between member nodes is allowed.  Each member can have as many of these
%cases superimposed on itself as necessary.  Just keep entering member
%number and the required input.
nfefcases=1; %Enter the number of separate fixed end force cases
mfef(1)=1;  %This says member one will have a some fixed end forces to consider
%mfef(2)=1;  %This says member there is another fixed end force case for member 1.
%mfef(3)=2;   %This says member 2 has a fixed end force case to consider.
fefcase(1)=2;  %This says fefcase 1 is a point load (code is 1 for point load);
%fefcase(2)=2;  %This says fefcase 2 is a distributed load (code is 2 for dist. load);
%fefcase(3)=2;  %This says fefcase 3 is a distributed load (code is 2 for dist. load);
fefdata(1,1)=-0.200/12; fefdata(1,2)=0;  %This says that fefcase 1 has a -5 force at 0.5 of member
% length from node 1;
%fefdata(2,1)=-2; fefdata(2,2)=0; %This says that fefcase 2 has a -2 dist load(fefdata(2,2)=0 this case).
%fefdata(3,1)=-1; fefdata(3,2)=0; %This says that fefcase 3 has a -1 dist load(fefdata(3,2)=0 this case).

%Zero all of the global fixed end forces
for i=1:nnodes*3;
    Fef(i)=0.0;
end;    
%Loop over the number of fixed end forces and add the fixed end forces into
%Fef one at a time.
for i=1:nfefcases;
    %zero the fixed end forces in local coordinates
    for j=1:6;
        feforce(j)=0.0;
    end;        
        
    %Determine which load case we must add in
    %Check if it is the point load case.
    if fefcase(i) == 1;         
        %add in the shear terms(negative point loads are downwards in local
        %coordinates).
        feforce(2)=-fefdata(i,1)*(1-fefdata(i,2)); %left shear
        feforce(5)=-fefdata(i,1)*(fefdata(i,2)); %right shear
        %left moment
        feforce(3)=-fefdata(i,1)*((1-fefdata(i,2))*L(mfef(i)))^2*fefdata(i,2)/L(mfef(i));
        %right moment
        feforce(6)=+fefdata(i,1)*(fefdata(i,2)*L(mfef(i)))^2*(1-fefdata(i,2))/L(mfef(i));
                
    %Check if it is the distributed load case.
    elseif fefcase(i) == 2;
        %add in the shear terms(negative point loads are downwards in local
        %coordinates).
        feforce(2)=-fefdata(i,1)*L(mfef(i))/2; %left shear
        feforce(5)=-fefdata(i,1)*L(mfef(i))/2; %right shear
        %left moment
        feforce(3)=-fefdata(i,1)*L(mfef(i))^2/12;
        %right moment
        feforce(6)=+fefdata(i,1)*L(mfef(i))^2/12;
       
    end;
    %Transform feforce to global coordinates
        %create T, the transformation matrix
        m=mfef(i);
        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;
        feft=T'*feforce';
        
        %add the forces into the global fixed end force matrix
        Fef(mconn(mfef(i),1)*3-2)=Fef(mconn(mfef(i),1)*3-2)+feft(1); 
        Fef(mconn(mfef(i),2)*3-2)=Fef(mconn(mfef(i),2)*3-2)+feft(4);       
        Fef(mconn(mfef(i),1)*3-1)=Fef(mconn(mfef(i),1)*3-1)+feft(2); 
        Fef(mconn(mfef(i),2)*3-1)=Fef(mconn(mfef(i),2)*3-1)+feft(5); 
        Fef(mconn(mfef(i),1)*3)=Fef(mconn(mfef(i),1)*3)+feft(3);
        Fef(mconn(mfef(i),2)*3)=Fef(mconn(mfef(i),2)*3)+feft(6);
end; %end loop over the fixed end force cases    

%***INPUT ENDS HERE ***

%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;

%Subtract the global fixed end forces from the global nodal forces
F=F-Fef;

%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;
    F(sup(i))=0.0; %zero nodal forces where the displacements are zero.
    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+Fef'

%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);
    
    %zero the fixed end forces in local coordinates
    for j=1:6;
        feforce(j)=0.0;
    end;        
    %Calculate the fixed end forces in local coordinates
    %Loop over the number of fixed end forces and add the fixed end forces into
%feforce one at a time.
for i=1:nfefcases;
    
    if mfef(i)== m;    
    %Determine which load case we must add in
    %Check if it is the point load case.
    if fefcase(i) == 1;         
        %add in the shear terms(negative point loads are downwards in local
        %coordinates).
        feforce(2)=feforce(2)-fefdata(i,1)*(1-fefdata(i,2)); %left shear
        feforce(5)=feforce(5)-fefdata(i,1)*(fefdata(i,2)); %right shear
        %left moment
        feforce(3)=feforce(3)-fefdata(i,1)*((1-fefdata(i,2))*L(mfef(i)))^2*fefdata(i,2)/L(mfef(i));
        %right moment
        feforce(6)=feforce(6)+fefdata(i,1)*(fefdata(i,2)*L(mfef(i)))^2*(1-fefdata(i,2))/L(mfef(i));
                
    %Check if it is the distributed load case.
    elseif fefcase(i) == 2;
        %add in the shear terms(negative point loads are downwards in local
        %coordinates).
        feforce(2)=feforce(2)-fefdata(i,1)*L(mfef(i))/2; %left shear
        feforce(5)=feforce(5)-fefdata(i,1)*L(mfef(i))/2; %right shear
        %left moment
        feforce(3)=feforce(3)-fefdata(i,1)*L(mfef(i))^2/12;
        %right moment
        feforce(6)=feforce(6)+fefdata(i,1)*L(mfef(i))^2/12;
       
    end;
end; %end checking if the fixed end force case is for the present member
end; %end loop over the fixed end force cases    
    
    %calculate the local element end forces
    display ('Local End Forces for element number:')
    m
    f=k*T*u'+feforce'
        
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=20000; %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


