%Code for implementing the example problem using MATLAB
% Ideas of assembly along with distretization will be implemented; 

clc; clear all;

%-------Input-----------------%
E = input('Youngs modulus of the material (in GPa); ');
L = input('Height of the tapered cyclinder (in m): ');
r0 = input('Radius of the tapered cyclinder at the fixed end (in m): ');
rL =input('Radius of the tapered cyclinder at the hanging end (in m): ');
n = input('Number of elements: ');
F = input('Force pulling down the body (in N): '); %Positive value of force will be treated as if it it pulling down

%-------Calculation---------%
Le = L/n;      %length of one element

x = zeros(n+1, 1);            %to store x-coordinates
A = zeros(n,1);               %to store area values
sigma = zeros(n,1);           %to store elemental stress values

k_e = zeros(2,2);             %Local 2 x 2 stiffness matrix
K_G = zeros(n+1,n+1);         %Global stiffness matrix of size DOF x DOF
F_G = zeros(n+1, 1);          %Global force matrix


for i = 1:n
    x_top = (i-1)*Le;                 % cordinate of top of the element
    x_bottom = i*Le;                  % cordinate of bottom of the element
    x(i+1) = x_bottom;                   %x for node 2
    r_top = r0-(r0-rL)/L*x_top;       %radius at the top of the element
    r_bottom = r0-(r0-rL)/L*x_bottom; %radius at the bottom of the element
    
    A(i) = pi*(r_top^2+r_bottom^2)/2;  %area of the element
    k(i) = (E*10^9)*A(i)/Le;           %spring constant of the bar element
    
    k_e = k(i)*[1 -1;-1 1];            %local stiffness matrix
    
    %------Assembly to global matrix----------%
   K_e = zeros(n+1, n+1);        %Local stiffness matrix with additional rows and columns of zeros
   K_e([i:i+1],[i:i+1])=k_e;
   K_G = K_G+K_e; 
end

%------------------Solution-----------------------%

F_G(n+1) = F;                     %Boundary condition

K_G_act = K_G([2:n+1],[2:n+1]);    %Since node 1 displacement is zero
F_G_act = F_G([2:n+1]);            %First equation will be the part of the constrant equation
U_act = inv(K_G_act)*F_G_act;      %Calculation of unknown displacements
U_tot = [0; U_act];                %all displacements
F_reaction = K_G(1,[1:n+1])*U_tot; %reaction force


% solving it analytically
xa = linspace(0,x_bottom, 1000);

for i = 1:1000
    U_a(i) = F*L^2/(E*10^9*pi()*(r0-rL))*(1/(L*r0-(r0-rL)*xa(i))-1/(L*r0));
end

plot(x,U_tot,'o'); %Plotting the FEM solution in open circular symbol
hold on;           %Let the graph be retained
plot(xa,U_a);      %Plotting the analytical solution in line
xlabel('Distance from fixed edge (m)');
ylabel('Displacement (m)');

