% One dimensional finite element program
% Frank Fisher, Stevens Institute of Technology
% Mod/Sim 435, Spring 2013
% Edited the generic 1D FEM code to do a specific Stepped-rod problem

clear; clc;
disp('Franks`s code for a specific problem');
disp('From the Figure, node 1 at the top, node 3 at the bottom'); disp('');

numele=2; 
numnod=numele+1;

x=[0, 25, 40];                  % x-coordinates of nodes [in]
node= [1:numele; 2:numele+1];   % node stores the nodes of all elements
area=(pi/4).*[1.25^2, 2^2];     % element area [in^2]
young=[435*10^3*ones(1,numele)]; % modulus of polystyrene = 435 ksi = 3 GPa

% support conditions, ifix(i)=1 if node i is fixed, else zero
ifix=[1,zeros(1,numele)];
ifix(numnod)=1;			%this will chance the last BC (fixed/free); 1 = fixed, 0 = free

% applied forces
force=[0,12*10^3,0];            % applied load [lb]. Note positive down
%
bigk=[zeros(numnod,numnod)];    % zero bigk matrix to prepare for assembly
%
% loop over elements
for e=1:numele
    length=x(node(2,e))-x(node(1,e));      % compute element length
    c=young(e)*area(e)/length;
    ke=[c,-c;-c,c];                        % compute element stiffness    
    % assemble ke into bigk
    bigk(node(1,e),node(1,e))=bigk(node(1,e),node(1,e))+ke(1,1);
    bigk(node(1,e),node(2,e))=bigk(node(1,e),node(2,e))+ke(1,2);
    bigk(node(2,e),node(1,e))=bigk(node(2,e),node(1,e))+ke(2,1);
    bigk(node(2,e),node(2,e))=bigk(node(2,e),node(2,e))+ke(2,2);
end

% support conditions (boundary conditions) - normal code
for n=1:numnod
    if (ifix(n) == 1)
        bigk(n,n)=1E+30;      force(n)=0;
    end
end
%
displace=force/bigk;        % solve stiffness equations

% compute stresses based on the nodel displacements
for e=1:numele
    length=x(node(2,e))-x(node(1,e));       % compute element length
    elong=displace(node(2,e))-displace(node(1,e));
    stress(2*e-1)=young(e)*elong/length;
    stress(2*e)=stress(2*e-1);
    xx(2*e-1)=x(node(1,e));
    xx(2*e)=x(node(2,e));
end

%break

% plot displacements and stresses
subplot(211), plot(x,displace,'*'); xlabel('position'), ylabel('displacement');
title('Displacement verses position along stepped rod');
subplot(212)
plot(xx,stress,'*-');; xlabel('position'), ylabel('stress');
title('Stress verses position along stepped rod');

disp('The displacements and stresses are ');
disp(displace);
disp(stress);

