%%%
%   
%    MATLAB FEM PROGRAM
%    By Anurag Gupta
%    Description:
%
%     Bar Under Axial Load
%     The external load is distributed along the length of 
%     the bar in  quadratic varying fashion: q(x) = cx^2,
%     where c is a constant. We have taken P=1. 

clear all;

%
% Define constants: 
% 
%  A: the cross section of the bar
%  E: Young's modulus
%
AE=1;
c=1;
P=1;  %Applied point load
Ln=1; %Length of the bar
u0=0; %dispacement boundary condition. u(0)= 0;
NumElem=20;  %number of elements
NumNode=NumElem + 1; %number of nodes

x=zeros(NumNode,1); %array for position of nodes
h=Ln/NumElem; %size of one element. All elements are taken to be of equal size

%Storing Node positions
tj=1;
for ti=0:h:Ln,
    x(tj)=ti;
    tj=tj+1;
end        

K=zeros(NumNode,NumNode); %Global Stiffness Matrix
R=zeros(NumNode,1); %Global Load vector
un=zeros(NumNode,1); %Computed Nodal Displacements

%running a loop for number of elements

for i = 1:NumElem,
    
     % We take i'th Element at a time 
     % And get the nodal positons for it
     xi_1=x(i);
     xi_2=x(i+1);
     
     % calculating Ke, ie local stiffness matrix
     ck = AE/h; 
     Ke = [ck -ck; -ck ck];
     
     %calculating Re, ie local load vector
     Re=zeros(2,1);
     
     Re(1,1)=(c/h)*(((xi_2^4)/12)-(xi_2*(xi_1^3)/3)+((xi_1^4)/4));
     Re(2,1)=(c/h)*(((xi_1^4)/12)-(xi_1*(xi_2^3)/3)+((xi_2^4)/4));
     % Assembling K, ie into global stiffness matrix
     
     K(i,i)=K(i,i)+Ke(1,1);
     K(i,i+1)=K(i,i+1)+Ke(1,2);
     K(i+1,i)=K(i+1,i)+Ke(2,1);
     K(i+1,i+1)=K(i+1,i+1)+Ke(2,2);
     
     % Assembling R, ie into global load vector
     
     R(i)=R(i)+Re(1);
     R(i+1)=R(i+1)+Re(2);
     
 end
 R(NumNode)=R(NumNode)+P;
 
 un(1) = u0; %enforcing boundar condition, at x=0
 un(2:NumNode)=K(2:NumNode,2:NumNode)\R(2:NumNode); %and determning rest of the nodal displacements
 
 %plotting
 
 figure;
x2=0:h:Ln;
 plot (x2,un,'-k')
 hold on; 
%plotting exact solution
x1=0:0.01:Ln;
u = (c./(AE)).*((-(c.*x1.^4)/12) + (P+((c.*Ln.^3)/3)).*x1);
plot(x1,u,'-r')

s = sprintf('FEM Analysis of Bar with an Axial Loading. Total number of Elements: %d', NumElem);
title(s)
title(s)
xlabel('u')
ylabel('x')
legend('FEM solution','Exact Solution',0)

%%%
%%%  End
%%%