%Mason Averill
%ME-480 Fall 2020, 11/20
%FEM Method Solution




%Input Parameters

l=1; %length of the beam in meters
t=0.05; %width of the beam in meters
e=200*10^9; %Young's Modulus of the steel utilized in Pa
density=7900; %density of the steel utilized in kg/m^3
gravity=9.81; %acceleration due to gravity
p=5000; %Point load applied at end of beam in N
number_of_elements_max=100; %specify number of elements to consider


number_of_elements=1;
FEM_end_deflection_store=zeros(1,number_of_elements_max);
FEM_end_slope_store=zeros(1,number_of_elements_max);
percent_error_deflection_store=zeros(1,number_of_elements_max);
while(number_of_elements<=number_of_elements_max)
%Discretize x
delta_x=l/number_of_elements;
x=0:delta_x:l;

%Optimum h(x) determined by optimizer: 0.063427138414689221796410276261629*x1 + 199/10000
%Now find h(x) from the wall towards the end of the beam

h_x=-0.063427138414689221796410276261629*x+199/10000+0.063427138414689221796410276261629;

%First lets find the area at each node

[rows,columns]=size(h_x);


A_x=zeros(rows,columns-1); %this will hold the area for each element
A_i=0;
A_j=0;
i=1;

while(i<columns)
 A_i=t*h_x(1,i);
 A_j=t*h_x(1,i+1);
 
 A_x(1,i)=(A_i+A_j)/2;   
    
 i=i+1;   
    
end 

%Now lets find I for each element

I_x=zeros(rows,columns-1);

i=1;

while(i<columns)
 I_x(1,i)=(t*((h_x(1,i)+h_x(1,i+1))/2)^3)/12;  
       
 i=i+1;   
end 


%Now lets assemble the keq for each element;


k_cell={};%create a cell array of all the k local matrices
k_local=zeros(6+3*(number_of_elements-1),6+3*(number_of_elements-1));%initialize the size of k local, in this case
%we are setting the size of k local equal to that of k global to make
%assembly of k global much easier at the end 


i=1;%used to index k local to be generated
A=0;
I=0;
start_row=0;
start_column=0;
l=delta_x;
while(i<=number_of_elements)
    start_row=(3*i)-2;
    start_column=start_row;
    
    A=A_x(1,i);%finds the area corresponding to the klocal in question
    I=I_x(1,i);%finds the second area moment of inertia for the klocal in question
    
    %now lets populate the k_local
    
    %populate first row
    
    k_local(start_row,start_column)=A*e/l;
    k_local(start_row,start_column+3)=-A*e/l;
    
    %populate second row
    k_local(start_row+1,start_column+1)=12*e*I/l^3;
    k_local(start_row+1,start_column+2)=6*e*I/l^2;
    k_local(start_row+1,start_column+4)=-12*e*I/l^3;
    k_local(start_row+1,start_column+5)=6*e*I/l^2;
    
    
    %populate third row
    k_local(start_row+2,start_column+1)=6*e*I/l^2;
    k_local(start_row+2,start_column+2)=4*e*I/l;
    k_local(start_row+2,start_column+4)=-6*e*I/l^2;
    k_local(start_row+2,start_column+5)=2*e*I/l;
    
    
    %populate fourth row
    k_local(start_row+3,start_column)=-A*e/l;
    k_local(start_row+3,start_column+3)=A*e/l;
    
    %populate fifth row
    k_local(start_row+4,start_column+1)=-12*e*I/l^3;
    k_local(start_row+4,start_column+2)=-6*e*I/l^2;
    k_local(start_row+4,start_column+4)=12*e*I/l^3;
    k_local(start_row+4,start_column+5)=-6*e*I/l^2;
    
    
    %populate sixth row
    k_local(start_row+5,start_column+1)=6*e*I/l^2;
    k_local(start_row+5,start_column+2)=2*e*I/l;
    k_local(start_row+5,start_column+4)=-6*e*I/l^2;
    k_local(start_row+5,start_column+5)=4*e*I/l;
    
    
    k_cell{1,i}=k_local;
 
    k_local=zeros(6+3*(number_of_elements-1),6+3*(number_of_elements-1));
  i=i+1;  
end 


%now lets assemble k Global

k_global=zeros(6+3*(number_of_elements-1),6+3*(number_of_elements-1));

i=1;

while(i<=number_of_elements)
   k_global=k_global+k_cell{1,i}; 
   
   i=i+1; 
end 

%k_global is now assembled

k_global_solve=k_global;
%now lets apply bc's and simplify 
k_global_solve(1,:)=0;
k_global_solve(2,:)=0;
k_global_solve(3,:)=0;
k_global_solve(:,1)=0;
k_global_solve(:,2)=0;
k_global_solve(:,3)=0;
k_global_solve(1,1)=1;
k_global_solve(2,2)=1;
k_global_solve(3,3)=1;

%now lets assemble the load matrix with b.c's applied

F=zeros((6+3*(number_of_elements-1)),1);

[rows, columns]=size(F);

F(rows-1,1)=-p;

u=k_global_solve^-1*F;

R=k_global*u-F;
l=1;

%Now lets find the analytical deflection
x1=0:delta_x:l;

h_x1=0.063427138414689221796410276261629*x1 + 199/10000;

analytical_deflection=-12*p/(e*t)*trapz(x1,x1.^2./h_x1.^3);

%Percent error between analytical and FEM method

percent_error=abs(analytical_deflection-u(rows-1,1))/abs(analytical_deflection)*100;

FEM_end_deflection_store(1,number_of_elements)=u(rows-1,1);
FEM_end_slope_store(1,number_of_elements)=u(rows,1);
percent_error_deflection_store(1,number_of_elements)=percent_error;
number_of_elements=number_of_elements+1;
end

%grab the deflection for each delta x

[rows,columns]=size(u);
[rowsx,columnsx]=size(x);
i=2;
j=1;
u_x=zeros(1,columnsx);
while(i<=rows)
 u_x(1,j)=u(i,1);   
    
 j=j+1;   
 i=i+3;   
end 

%grab the slope for each delta x


i=3;
j=1;
theta_x=zeros(1,columnsx);
while(i<=rows)
 theta_x(1,j)=u(i,1);   
    
 j=j+1;   
 i=i+3;   
end 



%show ideal h(x) function
figure (1)
plot(x,h_x);
xlabel('x (m)')
ylabel('h(x) (m)')
title('Optimum Height of the Beam from the Wall Towards the Free End')
grid on

%show that the beam weight is equal to the maximum allowable 
beam_weight=num2str(trapz(x,h_x)*density*gravity*t);

string_for_print=strcat('The weight of the beam is ',{' '},beam_weight,' Newtons');

disp(string_for_print);

elements_for_plot=1:1:number_of_elements_max;


figure (2)

plot(elements_for_plot,FEM_end_deflection_store)
title('Deflection at the End of the Beam Determined by FEM method vs Number of Elements Considered')
xlabel('Number of elements')
ylabel('Deflection at the End of the Beam')
grid on

FEM_end_slope_store=FEM_end_slope_store*180/pi();

figure (3)
plot(elements_for_plot,FEM_end_slope_store)
xlim([2,number_of_elements_max])
title('Slope at the Free End of the Beam Determined by FEM method vs Number of Elements Considered')
xlabel('Number of Elements')
ylabel('Slope at the Free End of the Beam (degrees)')
grid on

figure (4)
plot(elements_for_plot,percent_error_deflection_store)
title('Percent Error of FEM Method Compared to Analytical Solution vs Number of Elements Considered')
xlabel('Number of Elements')
ylabel('Percent error (%)')
xlim([10,100])
grid on

figure (5)
plot(x,u_x)
title('Deflection of Beam From the Wall Towards the Free End')
xlabel('x (m)')
ylabel('Deflection of Beam (m)')
grid on

figure (6)
plot(x,theta_x)
title('Slope of Beam vs Position in x')
xlabel('x (m)')
ylabel('Slope of Beam (rad)')
grid on

theta_x=theta_x*180/pi();

figure (7)
plot(x,theta_x)
title('Slope of Beam From the Wall Towards the Free End')
xlabel('x (m)')
ylabel('Slope of Beam (degrees)')
grid on