function [] = beam1D
%------------------------------------------------------------
% Beam elements with 2 Degrees Of Freedom at each node.
% Imad Barsoum & Jonas Faleskog /KTH Solid Mechanics/2006-Feb-8/
%----------------------------------------------------------------
% 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
%      element: Columns: 1 = 1st node #, 2 = 2nd node #, 3 = E modulus,
%                        4 = Area, 5 = momen of intertia Iz
%      D:   global displacement vector
%      F:   global force vector
%------------------------------------------------------------
% 1. Read in-put data from file = beam1D.inp
%       The global force vector is established here!
fid = fopen('beam1D.inp','r');
[node,element,Dzero,Fs,elem_force,Nn,Ne] = read_input(fid)
fclose(fid);
%------------------------------------------------------------
% 2. Assembly of (i) global stiffness matrix of the system of trusses
%       and (ii) global force vector.
%       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
%           Ksystem(i,j)  dimension: No_of_Eqns x No_of_Eqns
%
No_of_Eqns = 2*Nn;
Fd = zeros(No_of_Eqns,1);
Ksystem = zeros(No_of_Eqns,No_of_Eqns);
for e = 1 : Ne
    Enode1 = round( element(e,1) ); Enode2 = round( element(e,2) ); 
    Eemod = element(e,3); EIz = element(e,4);
    [Ke] = ElemStiffMatrix(Eemod,EIz,node(Enode1),node(Enode2));
    [Fde] = ElemDistForces(elem_force(e,:),node(Enode1),node(Enode2));
    ekv(2) = 2*Enode1;
    ekv(1) = ekv(2) - 1;
    ekv(4) = 2*Enode2;
    ekv(3) = ekv(4) - 1;
    for i = 1 : 4
        k = ekv(i);
        Fd(k) = Fd(k) + Fde(i);
        for j = 1 : 4
            l = ekv(j);
            Ksystem(k,l) = Ksystem(k,l) + Ke(i,j);
        end
    end
end
F = Fs + Fd;
%------------------------------------------------------------
% 3. Enforce zero displacment/rotation 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 = Ksystem;
for i = 1 : No_of_Eqns
    if( Dzero(i) == 1 )
        Kmod(i,:) = 0;
        Kmod(:,i) = 0;
        Kmod(i,i) = 1.0;
        F(i) = 0;
    end
end
%------------------------------------------------------------
% 4. Solve the Equation system
D = Kmod\F
%------------------------------------------------------------
% 5. Compute the reaction forces/moments (this is simply done by multiplying
%    the displacement vector by the full system stiffness matrix)
F = Ksystem*D
%------------------------------------------------------------
% 6. Plot the deformed and the undeformed configuration
dmag = input('Give magnification factor for displacements in plot: ');
plot_result(node,element,Nn,Ne,D,F,dmag);
%------------------------------------------------------------
% 7. Generating the output file beam1D.out where the results are listed
gen_output(node,element,D,F,Fd,Dzero);
%*************************************************************
%           FUNCTIONS     and     SUBROUTINES
%*************************************************************
% Compute the element stiffness matrix 
function [Ke] = ElemStiffMatrix(E,Iz,node1,node2)
%  Retrieve coordinates of element nodes
x1 = node1(1);  x2 = node2(1);
%  Element length
le = (x2-x1); a = le/2;
%  Stiffness of the beam element
beam_stiffness = E*Iz/(2*a^3);
%  Element stiffness matrix
Ke = [ 3     3*a     -3     3*a   ; ...
       3*a   4*a^2   -3*a   2*a^2 ; ...
      -3    -3*a      3    -3*a   ; ...
       3*a   2*a^2   -3*a   4*a^2 ;];
Ke = beam_stiffness*Ke;
%*************************************************************
% Compute consistent nodal forces given a distributed force acting along
%  the element on the form: fd = q0 + q1*(x/le) + q2*(x/le)^2
function [Fd] = ElemDistForces(elem_force,node1,node2)
Fd = zeros(4,1);
if (elem_force(1) ~= 0)
% Retrieve coordinates of element nodes
    x1 = node1(1);  x2 = node2(1);
% Element length
    le = (x2-x1); a = le/2;
% Compute node forces
    qi = elem_force(2); qj = elem_force(3);
    Fi =  a*qi * [0.7 ; 0.2*a ; 0.3 ; -2*a/15];
    Fj =  a*qj * [0.3 ; 2*a/15 ; 0.7 ; -0.2*a];
    Fd =  Fi + Fj;
end
%*************************************************************
%  This routine generates a plot of deformed and undeformed configuration
function plot_result(node,element,Nn,Ne,D,F,dmag)
% Undeformed configuration
%  Plot nodes
for i = 1 : Nn
    x = node(i);  y = 0.0; subplot(3,1,1); plot(x,y,'ob'); hold on
end
%  Plot elements
for i = 1 : Ne
    Enode1 = round(element(i,1));  Enode2 = round(element(i,2)); 
    x = [node(Enode1,1),node(Enode2,1)];   y = [0,0];   line(x,y);
end
% Deformed configuration
%  Plot nodes
xi = linspace(-1,1,100);
w=[]; x_coord=[]; M=[]; T=[];
for i = 1 : Ne
    Emod = element(i,3);  Iz = element(i,4);
    Enode1 = round(element(i,1));  Enode2 = round(element(i,2)); 
    x1 = node(Enode1,1);   x2 = node(Enode2,1);  le = abs(x2-x1); a = le/2;
    x_e = linspace(x1,x2,length(xi));  de = [D(2*Enode1-1:2*Enode1);D(2*Enode2-1:2*Enode2)];
    for j = 1 : length(xi)
        % Determine the Displacement
        N1 = (2 - 3*xi(j) + xi(j)^3)./4;
        N2 = (1 - xi(j) - xi(j)^2 + xi(j)^3)*a./4;
        N3 = (2 + 3*xi(j) - xi(j)^3)./4;
        N4 = (-1 - xi(j) + xi(j)^2 + xi(j)^3)*a./4;
        N = [N1 N2 N3 N4];   w_xi(j) = N*de;
        % Determine the Moments
        N1_2 =  3*xi(j)./2;
        N2_2 =  a*(3*xi(j)-1)./2;
        N3_2 = -3*xi(j)./2;
        N4_2 =  a*(1+3*xi(j))./2;
        N_biss = [N1_2 N2_2 N3_2 N4_2];  M_xi(j) = -Emod*Iz*N_biss*de/a^2;
        % Determine the Forces (shear forces)
        N_triss = 3*[1 a -1 a]./2;  T_xi(j) = -Emod*Iz*N_triss*de/a^3;
    end
    x_coord = [x_coord x_e]; w = [w dmag*w_xi];  M = [M dmag*M_xi];
    T = [T dmag*T_xi];
end
subplot(3,1,1);
plot(x_coord,w); ylabel('w'); hold on;
subplot(3,1,2); 
plot(x_coord,M); ylabel('M'); hold on;
subplot(3,1,3);
plot(x_coord,T); ylabel('T'); hold on;
%axis([min(x_coord) max(x_coord) min(v_disp) max(v_disp)]);
%*************************************************************
% This routine reads in-put data from the file "beam1D.inp"
function [node,element,Dzero,Fs,elem_force,Nn,Ne] = read_input(fid)
data = [];
% Node numbers and coordinates
data = str2num(fgetl(fid));
data = str2num(fgetl(fid));
while (length(data) > 0)
    Nnum = round(data(1));
    node(Nnum,1) = data(2);
    data = str2num(fgetl(fid));
end
Nn = length(node(:,1));
%
% Elements (connectivity, E & A)
data = str2num(fgetl(fid));
while (length(data) > 0)
    Enum = round(data(1));
    element( Enum, 1:4) = data(2:5);
    data = str2num(fgetl(fid));
end
Ne = length(element(:,1));
%
% Zero disp./rotation Boundary Conditions
Dzero = zeros(2*Nn,1);
data = str2num(fgetl(fid));
while (length(data) > 0)
    if (data(2) == 1)
        Dzero(2*data(1)-1) = 1;
    elseif (data(2) == 2)
        Dzero(2*data(1)) = 1;
    end
    data = str2num(fgetl(fid));
end
%
% Applied point forces/moments
Fs = zeros(2*Nn,1);
data = str2num(fgetl(fid));
while (length(data) > 0)
    if (data(2) == 1)
        Fs(2*data(1)-1)  = data(3);
    elseif (data(2) == 2)
        Fs(2*data(1)) = data(3);
    end
    data = str2num(fgetl(fid));
end
%
% Applied distributed force (perpendicular to the beam-element)
elem_force = zeros(Ne,3);
data = str2num(fgetl(fid));
while (length(data) > 0)
    elem_force(data(1),:) = data;
    data = str2num(fgetl(fid));
end
return
%*************************************************************
% This routine generates an output file where the results are listed
function gen_output(node,element,D,F,Fb,Dzero)
fid = fopen('beam1D.out','w+');
t = clock;
fprintf(fid,'%s\n','****************************************************');
fprintf(fid,'%s%4.2f\n',['File: beam1D.out   Date: ',date,'   Time: '],t(4) + .01*t(5));
fprintf(fid,'%s\n','****************************************************');
fprintf(fid,'%s\n','------------- External Load Vectors ----------------');
fprintf(fid,'%s\n','Node           Fz                  My');
for i = 1 : length(node(:,1))
    if ( Dzero(2*i-1)==0 & Dzero(2*i)==0 ) 
       fprintf(fid,'%u%s%e%s%e\n',i,'    ',F(2*i-1),'    ',F(2*i));
    elseif  ( Dzero(2*i-1)==0 & Dzero(2*i)==1 )
       fprintf(fid,'%u%s%e%s\n',i,'    ',F(2*i-1),'              - ');
    elseif  ( Dzero(2*i-1)==1 & Dzero(2*i)==0 )
       fprintf(fid,'%u%s%e\n',i,'              -         ',F(2*i));
    elseif  ( Dzero(2*i-1)==1 & Dzero(2*i)==1 )
       fprintf(fid,'%u%s\n',i,'              -                   -');
    end
end
fprintf(fid,'%s\n','------------- Nodal DOF ------------------');
fprintf(fid,'%s\n','Node           Dz                  Thetay');
for i = 1 : length(node(:,1))
    fprintf(fid,'%u%s%e%s%e\n',i,'    ',D(2*i-1),'    ',D(2*i));
end
fprintf(fid,'%s\n','------------- Reaction Forces/Moments ----');
fprintf(fid,'%s\n','Node           Rz                  My');
R=F-Fb;
for i = 1 : length(node(:,1))
    if ( Dzero(2*i-1)==1 & Dzero(2*i)==1 ) 
       fprintf(fid,'%u%s%e%s%e\n',i,'    ',R(2*i-1),'    ',R(2*i));
    elseif  ( Dzero(2*i-1)==1 & Dzero(2*i)==0 )
       fprintf(fid,'%u%s%e%s\n',i,'    ',R(2*i-1),'              - ');
    elseif  ( Dzero(2*i-1)==0 & Dzero(2*i)==1 )
       fprintf(fid,'%u%s%e\n',i,'              -         ',R(2*i));
    elseif  ( Dzero(2*i-1)==0 & Dzero(2*i)==0 )
       fprintf(fid,'%u%s\n',i,'              -                   -');
    end
end
fclose(fid);
disp(' *** The result is listed in the output file beam1D.out ***');
%*************************************************************