% One dimensional finite element program
% Frank Fisher, Stevens Institute of Technology
% Simulation and Modeling 345
%
% INPUT DATA
% this data automatically describes a line of 2-node elements
% with support conditions at both ends and a uniform load; the
% number of elements can be varied by changing numele

clear;
clc;
close all;                      %closes all open Figure windows
disp('This is a general Matlab code for a 1D rod.');
disp('This program can be edited to solve a number of suitable problems');
numele=2;                       % number of elements
numnod=numele+1;                % number of nodes (special case for 1D problem)
x= 0:1/numele:1;                % x-coordinates of nodes
node= [1:numele; 2:numele+1];   % node stores the nodes of all elements
area=[1*ones(1,numele)];        % area of each element
young=[1E7*ones(1,numele)];     % Young's modulus of each element
% support conditions, ifix(i)=1 if node i is fixed, else zero
ifix=[1,zeros(1,numele)];
ifix(numnod)=0;			        %change the BC on last node; 1 = fixed, 0 = free		
force=[1/numele*ones(1,numnod)];% applied external forces on nodes
%
% BELOW IS THE MACHINERY OF THE FEM METHOD
% DO NOT CHANGE - SKIP TO THE PLOTTING PART
bigk=[zeros(numnod,numnod)];    % zero bigk matrix to prepare for assembly
%
for e=1:numele                          % loop over elements
    length=x(node(2,e))-x(node(1,e));   % compute element length
    c=young(e)*area(e)/length;          % "spring stiffness" for each rod
    ke=[c,-c;-c,c];                     % compute element stiffness%
    % now 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)
for n=1:numnod
    if (ifix(n) == 1)
        bigk(n,n)=1E+30;                % math trick to account for fixed BC
        force(n)=0;
    end
end
%
disp=force/bigk;                        % SOLVE [F]=[K][D] FOR THE NODAL DISPLACEMENTS
%
subplot(211), plot(x,disp,'*')          % plot displacements
title('Displacements calculated via FEM');
ylabel('displacement')

% compute stresses

for e=1:numele                              % this loops over the number of elements
    length=x(node(2,e))-x(node(1,e));       % compute element length
    elong=disp(node(2,e))-disp(node(1,e));  % elongation of each element
    stress(2*e-1)=young(e)*elong/length;    % to plot, calculate the "stress" at each node
    stress(2*e)=stress(2*e-1);              % stress in each element is uniform
    xx(2*e-1)=x(node(1,e));
    xx(2*e)=x(node(2,e));
end

subplot(212)
plot(xx,stress,'*-')  % plot the FEM stresses (constant stress element)

title('Stresses calculated via FEM');
xlabel('x position')
ylabel('Stresses')
