%#########################################################################
%%%%% Version 2.0 #######  MOHAMED ABDOU MAHRAN KASEM ##################
% Ph.D. Aerospace Engineering Department, Cairo Uiniversity
% Egypt (17th Oct. 2015)
function [K, Vel] =truss_stiff(elem,coord,E,A, nel, edof, sdof)
%--------
% purpose
%--------
%  The Global Stiffness Matrix
%----------------------
% Variables Description
%----------------------
%-----------------------------------------------------------------------
%  Variable      |                     Description
%-----------------------------------------------------------------------
%   coord        |                  nodes coordinates matrix
%    elem        |                 element conectivity matrix
%     nel        |                 number of elements
%     E          |                  Elastic Modulus
%     nu         |                   Poison ratio
%     k          |  element stiffness matrices its a 3-Dim matrix 
%     F          |                 Global force vector
%     D          |               Global displacement vector
%     B          |             the strain-displacement matrix
%     sdof       |             total system dof
%     edof       |            total element dof 
%     eldof      |            element dof as a vector
%     k          |     element stiffnes matrix in global coordinates
%-----------------------------------------------------------------------

%---------------------------------------
% Initialization of matrices and vectors
%---------------------------------------
    k = zeros(edof,edof);
    K = zeros(sdof, sdof);
    Vel = zeros(1,nel);
for i = 1:nel
    % element nodes
    ncon(1) = elem(i,1);  % 1st connected node
    ncon(2) = elem(i,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
    Vel(i)  = Le*A(i);              % element volume
    theta   = atan(dy/dx);
    
    l = cos(theta); m = sin(theta); %direction cosines
    C = (E*A(i))/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);
      
     % Assembly of the Global Stiffness matrix
      eldof          = [2*elem(i,:)-1; 2*elem(i,:)];
      K(eldof,eldof) = K(eldof,eldof) + k;
end