% 2D TRUSS PROGRAM
% Spefically written for the sample problem for the 2D truss
%
clear;                                  %clear previous runs

%**********INPUT DATA***************************************

numele =2;                              %number of elements
numnod =3;                              %number of nodes
x= [0,0,-10];                   %nodal x-coordinates
y= [0,10,0];                   %nodal y-coordinates

%%%connectivity of elements to nodes
%element I is connected to nodes 1&2
%element 2 is connected to nodes 2&3
%element 3 is connected to nodes 3&4
%etc
node= [1,3;
        2,2];

%Note: the size of 'young' and 'area' must match numele
young=1e6*[1 1];        % Young's modulus
area=1*[1 1];        % cross-sectional areas

%Nodal Constraints
ifix=[zeros(2,numnod)];         %initialize ifix
%ifix=0 is a free DOF, ifix=1 is a constrained DOF
%in this example nodes 1 and 6 fully constrained; others free

ifix=[1,0,1;
       1,0,1];

%Applied Force Vector
% force=[Fx1;Fy1;Fx2;Fy2;Fx3;Fy3;Fx4;Fy4;Fx5;Fy5;Fx6;Fy6]
% input only where known, leave as zero in all other slots
force=[0;0;-500;380;0;0];

%       END OF INPUT DATA
% ********************************************************

% zero bigk matrix to prepare for assembly
ndof=2;
numeqns=numnod*ndof;
bigk=[zeros(numeqns)];

%  loop over elements

for e=1:numele          % compute element stiffness
     nlink=2;            % nlink=number of nodes per element
     x21=x(node(2,e))-x(node(1,e));
          y21=y(node(2,e))-y(node(1,e));
     length=sqrt(x21^2+y21^2);
     length_print(e)=length;        % frank added to troubleshoot
     c=x21/length;
     s=y21/length;
     ke=area(e)*young(e)/length*[c^2, c*s,-c^2,-c*s;
                           c*s, s^2,-c*s,-s^2;
                          -c^2,-c*s, c^2, c*s;
                          -c*s,-s^2, c*s, s^2];

         % assemble ke into bigk

         n1=ndof-1;
    for i=1:nlink;
      for j=1:nlink;
        rbk=ndof*(node(i,e)-1)+1;
        cbk=ndof*(node(j,e)-1)+1;
        re=ndof*(i-1)+1;
        ce=ndof*(j-1)+1;

bigk(rbk:rbk+n1,cbk:cbk+n1)=bigk(rbk:rbk+n1,cbk:cbk+n1)+ke(re:re+n1,ce:ce+n1);
      end                %j loop
         end             %i loop
end                     %e loop

disp('Unconstrained Stiffnes Matrix is:')
disp(bigk)
disp('Hit any key to continue')
pause;

% impose support conditions (boundary conditions)
for n=1:numnod
      for j=1:ndof
          if (ifix(j,n) == 1)
            m=ndof*(n-1)+j;
            dummy=bigk(m,m);
            bigk(m,:)=zeros(1,numeqns);
            bigk(:,m)=zeros(numeqns,1);
            bigk(m,m)=1;
                    force(m)=0;  %zero appropriate rows of force
          end
     end
  end

% solve stiffness equations
displace=bigk\force                 %nodal displacements

% compute stresses
for e=1:numele
     n1=2*node(1,e)-1;           %get nodes of the element
     n2=2*node(2,e)-1;
     de(1:2)=displace(n1:n1+1);      %nodal disp for first node
     de(3:4)=displace(n2:n2+1);      %nodal disp for second node
     x21=x(node(2,e))-x(node(1,e));
     y21=y(node(2,e))-y(node(1,e));
     lengthvec(e)=sqrt(x21^2+y21^2);
     length=lengthvec(e);
     c=x21/length;
     s=y21/length;
     S=[-c,-s,c,s];
     stress(e)=(young(e)/length)*S*de';
end

stresscol=stress'               %print stresses in elements
