% One dimensional finite element program
% Frank Fisher, Stevens Institute of Technology
% Mod/Sim 435, Spring 2013
% Edited the standard 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=3; 
numnod=numele+1;

x=[0, 5, 30, 50]*(1/100);       % x-coordinates of nodes [cm]
node= [1:numele; 2:numele+1];   % node stores the nodes of all elements
area= (1/100)^2*[1, 0.5, 0.25]; % element area [m^2];convert from cm
young=10^9*[0.3, 6, 12.8];      % moduli in Pa

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

force=[0,50, 100, 80];            % applied forces [in N]. 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


for n=1:numnod;    % enforce support conditions (boundary conditions)
    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

% 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);

