%#########################################################################
%%%%% Version 2.0 #######  MOHAMED ABDOU MAHRAN KASEM ##################
% Ph.D. Aerospace Engineering Department, Cairo Uiniversity
% Egypt (17th Oct. 2015)
function [del, fel_local, sigma]=dis_truss(D, E,A, nel, elem, edof, coord)

% Purpose
% to calculate the element displacement and reaction forces

%----------------------
% Variables Description
%----------------------
%-----------------------------------------------------------------------
%  Variable      |                     Description
%-----------------------------------------------------------------------
%    elem        |                 element conectivity matrix
%     nel        |                 number of elements
%     F          |                 Global force vector
%     D          |               Global displacement vector
%     del        |            the element nodal displacement matrix
%     fel        |            the element nodal reaction matrix
%    dof         |                    all dofs
%    edof        |                  total element dof
%    eldof       |                   elemet dof vector 
%    theta       |                element orientation angle
%-----------------------------------------------------------------------
%---------------------------------------
% Initialization of matrices and vectors
%---------------------------------------
del       = zeros(edof, nel);
fel       = zeros(edof, nel);
fel_local = zeros(2, nel);
sigma     = zeros(1, nel);
for ie = 1:nel
    % element nodes
    ncon(1) = elem(ie,1);  % 1st connected node
    ncon(2) = elem(ie,2);  % 2nd connected node
    % coordinates element's nodes
    x1 = coord(ncon(1),1); y1 = coord(ncon(1),2);
    x2 = coord(ncon(2),1); y2 = coord(ncon(2),2);
    % element length and orientation
    dx      = x2-x1;
    dy      = y2-y1;
    Le      = sqrt(dx^2+dy^2);      % element length
    theta   = atan(dy/dx);
    
    l = cos(theta); m = sin(theta); %direction cosines
    T = [l m 0 0; 0 0 l m];
    % Element stiffness matrix
    C = (E*A(ie))/Le;
    
    % Element stiffness matrix in global coordinates
      k(1,1)=l^2*C;
      k(2,1)=l*m*C;    k(2,2)=m^2*C;
      k(3,1)=-l^2*C;   k(3,2)=-l*m*C;    k(3,3)=l^2*C;
      k(4,1)=-l*m*C;   k(4,2)=-m^2*C;    k(4,3)=l*m*C;    k(4,4)=m^2*C;
      
                       k(1,2)=k(2,1);    k(1,3)=k(3,1);   k(1,4)=k(4,1);
                                         k(2,3)=k(3,2);   k(2,4)=k(4,2);
                                                          k(3,4)=k(4,3);
    
    eldof     = [2*elem(ie,1)-1;2*elem(ie,1);2*elem(ie,2)-1;2*elem(ie,2)];
    del(:,ie) = D(eldof);
    fel(:,ie) = k*del(:,ie);
    fel_local(:,ie) = T*fel(:,ie);
    sigma(ie) = fel_local(2,ie)/A(ie); 
end