function [f,u] = spring2D
%------------------------------------------------------------
% 2D spring elements with 2 Degrees Of Freedom at each node.
% Imad Barsoum & Jonas Faleskog /KTH Solid Mechanics/2005-Jan-27/
%----------------------------------------------------------------
% Variables:
%    Problem size
%       Nn:         total number of nodes
%       Ne:         total number of elements
%       No_of_Eqns: total number of equations (degrees of freedom)
%    Vectors & Matrices:
%      node:    Columns: 1 = node #, 2 = x-coord, 3 = y-coord
%      element: Columns: 1 = element #,  2 = 1st node #, 
%                        3 = 2nd node #, 4 = spring constant
%------------------------------------------------------------
% 1. Read in-put data from file = spring2D.inp
%       The global force vector is established here!
fid = fopen('spring2D.inp','r');
[node,element,u,f,Nn,Ne] = read_input(fid);
fclose(fid);
%------------------------------------------------------------
% 2. Assembly of global stiffness matrix of the system of springs.
%       Loop over all elements (compute the element stiffness matrix
%       for each element in the global coordinate system and add it
%       to the system stiffness matrix).
%       Define the dimension of the matrix and set it to zero.
%
%       Element stiffnes matrix in global coordinate system
%           Ke(i,j)  dimension: 4x4 
%       System stiffnes matrix in global coordinate system
%           Ks(i,j)  dimension: No_of_Eqns x No_of_Eqns
%
No_of_Eqns = 2*Nn;
Ks = zeros(No_of_Eqns,No_of_Eqns);
for e = 1 : Ne
    [Ke] = ElemStiffMatrix(element(e,4),node(element(e,2),:),node(element(e,3),:));
    ekv(2) = 2*element(e,2);
    ekv(1) = ekv(2) - 1;
    ekv(4) = 2*element(e,3);
    ekv(3) = ekv(4) - 1;
    for i = 1 : 4
        k = ekv(i);
        for j = 1 : 4
            l = ekv(j);
            Ks(k,l) = Ks(k,l) + Ke(i,j);
        end
    end
end
%------------------------------------------------------------
% 3. Enforce zero displacment boundary conditions by setting the
%      appropriate columns and rows to zero in the system stiffness
%      matrix. Actual reduction of the system matrix can be avoided
%      by setting the diagonal term of constrained D.O.F. to unity.
Kmod = Ks;
for i = 1 : No_of_Eqns
    if( u(i) == 0 )
        Kmod(i,:) = 0;
        Kmod(:,i) = 0;
        Kmod(i,i) = 1.0;
    end
end
%------------------------------------------------------------
% 4. Solve the Equation system
u = Kmod\f
%------------------------------------------------------------
% 5. Compute the reaction forces (this is simply done by multiplying
%        the displacement vector by the full system stiffness matrix)
f = Ks*u
%------------------------------------------------------------
% 6. Plot the deformed and the undeformed configuration
dmag = input('Give magnification factor for displacements in plot: ')
plot_result(node,element,Nn,Ne,u,f,dmag);
%*************************************************************
%           FUNCTIONS     and     SUBROUTINES
%*************************************************************
% Compute the element stiffness matrix in the global coordinate
%    system for the element in question. 
function [Ke] = ElemStiffMatrix(spring_stiffness,node1,node2)
%  Retrieve coordinates of element nodes
x1 = node1(2);  x2 = node2(2);  y1 = node1(3);  y2 = node2(3);
%  Element length
le = sqrt( (x2-x1)^2 + (y2-y1)^2 );
%  Direction cosines
c = (x2-x1) / le;
s = (y2-y1) / le;
k_c2 = c*c * spring_stiffness;
k_s2 = s*s * spring_stiffness;
k_cs = c*s * spring_stiffness;
%  Element stiffness matrix transformed to global coordinate system
Ke = [ k_c2   k_cs   -k_c2   -k_cs ; ...
       k_cs   k_s2   -k_cs   -k_s2 ; ...
      -k_c2  -k_cs    k_c2    k_cs ; ...
      -k_cs  -k_s2    k_cs    k_s2 ;];
%*************************************************************
%  This routine generates a plot of deformed and undeformed configuration
function plot_result(node,element,Nn,Ne,u,f,dmag)
hold on
% axis(1.1*[min(node(:,2)) max(node(:,2)) min(node(:,3)) max(node(:,3))]);
% Undeformed configuration
%    Plot nodes
for i = 1 : Nn
    plot(node(i,2),node(i,3),'ob')
end
%    Plot elements
for i = 1 : Ne
    x = [node((element(i,2)),2),node((element(i,3)),2)];
    y = [node((element(i,2)),3),node((element(i,3)),3)];
    line(x,y)
end
% Deformed configuration
%    Plot nodes
for i = 1 : Nn
    plot(node(i,2)+dmag*u(2*i-1),node(i,3)+dmag*u(2*i),'or')
end
%    Plot elements
for i = 1 : Ne
    x = [node((element(i,2)),2) + dmag*u(2*element(i,2)-1),...
         node((element(i,3)),2) + dmag*u(2*element(i,3)-1)];
    y = [node((element(i,2)),3) + dmag*u(2*element(i,2)),...
         node((element(i,3)),3) + dmag*u(2*element(i,3))];
    line(x,y)
end
%*************************************************************
% This routine reads in-put data from the file "spring2D.inp"
function [node,element,u,f,Nn,Ne] = read_input(fid)
data = [];
% Node numbers and coordinates
data = str2num(fgetl(fid));
data = str2num(fgetl(fid));
while (length(data) > 0)
    node(data(1),:) = data;
    data = str2num(fgetl(fid));
end
Nn = length(node(:,1));
% Elements (connectivity)
data = str2num(fgetl(fid));
while (length(data) > 0)
    element(data(1),:) = data;
    data = str2num(fgetl(fid));
end
Ne = length(element(:,1));
% Zero displacement Boundary Conditions
u = ones(2*Nn,1);
data = str2num(fgetl(fid));
while (length(data) > 0)
    if (data(2) == 1)
        u(2*data(1)-1)  = 0.0;
    elseif (data(2) == 2)
        u(2*data(1)) = 0.0;
    end
    data = str2num(fgetl(fid));
end
% Applied point forces
f = zeros(2*Nn,1);
data = str2num(fgetl(fid));
while (length(data) > 0)
    if (data(2) == 1)
        f(2*data(1)-1)  = data(3);
    elseif (data(2) == 2)
        f(2*data(1)) = data(3);
    end
    data = str2num(fgetl(fid));
end
return
%*************************************************************