% Bar Under linearly distributed load and zero Axial Load
% The load is distributed along the length of the bar in linearly varying fashion: q = cx
% c is a constant. We have taken P=0. That is, there is no axial load at the end of the bar

clear all;

% constants
AE=1;
c=1;
Ln=1; %Length of the bar
u0=0; %dispacement boundary condition. u(0)= 0;
NumElem=2;  %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 = (c.*(h.^2)/6).*[ 3.*i - 2 ; 3.*i - 1];
     
     % 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
 
 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; 
%ploting exact solution
x1=0:0.01:Ln;
u = (c./(6.*AE)).*(3.*(Ln.^2).*x1 - (x1.^3));
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)