%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% finite element program for euler bernoulli beam elements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% example: 5 euler bernoulli beam elements, 6 nodes
%
%       1       2       3       4       5         elements
%   *-------*-------*-------*-------*-------*
%   1       2       3       4       5       6     nodes
% (1,2)   (3,4)   (5,6)   (7,8)   (9,10) (11,12)  degrees of freedom
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input data and connectivity matrix edof (el vs ndof)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all

disp('*** FINITE ELEMENT ANALYSIS OF EULER BERNOULLI BEAM ***');
EI     = input('bending stiffness EI = ');
length = input('beam length       l  = ');
loadp  = input('line load         p  = ');
loadS  = input('slear force       S  = ');
loadM  = input('bending moment    M  = ');

nel  = 5;
node = 6;

K  = zeros(2*node,2*node);
F  = zeros(2*node,1);
le = length/nel; 

edof = [ 1  1  2  3  4;
         2  3  4  5  6; 
         3  5  6  7  8; 
         4  7  8  9 10; 
         5  9 10 11 12];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% element stiffness matrix and element load vector 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ke = 2*EI/le^3*[    6   3*le    -6     3*le  ;
                 3*le   2*le^2  -3*le    le^2;
                   -6  -3*le     6    -3*le  ;
				 3*le     le^2  -3*le  2*le^2];

Fe  = le*loadp/2 * [ 1  ; le/6;  1  ;   -le/6];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% assembly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nie,n]=size(edof);
t=edof(:,2:n);
for e = 1:nie
	for i=1:(n-1)
		F(t(e,i))=F(t(e,i))+Fe(i);
        for j=1:(n-1)
            K(t(e,i),t(e,j)) = K(t(e,i),t(e,j))+Ke(i,j);
        end
	end    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% boundary conditons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bc = [1 0;  2 0];                   % left : uy(0) = 0, theta(0) = 0
F(2*node-1) = F(2*node-1) + loadS;  % right: S(L)
F(2*node)   = F(2*node)   - loadM;  % right: M(L)

%bc = [1 0; 2 0; 2*node-1 0; 2*node 0] % both: uy(0,L) = 0, theta(0,L) = 0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% partition system and solve
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nd,nd]=size(K);
fdof=[1:nd]';

d=zeros(size(fdof));
R=zeros(size(fdof));

pdof = bc(:,1); dp=bc(:,2);  fdof(pdof)=[];
s = K(fdof,fdof)\(F(fdof)-K(fdof,pdof)*dp);
d(pdof) = dp;  
d(fdof) = s; 
R = K*d-F; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% postprocessing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xi = -.96:0.04:1;
% cubic hermitian shape function
N10 =      (2-3*xi+xi.^3)/4;
N20 =      (2+3*xi-xi.^3)/4;
N11 = ( 1-xi-xi.^2+xi.^3)/4;
N21 = (-1-xi+xi.^2+xi.^3)/4;

% deflection elements 1,2,3,4,5
w(:,1) = N10'*d(1) + N11'*d( 2)*le/2 + N20'*d(3) + N21'*d( 4)*le/2;
w(:,2) = N10'*d(3) + N11'*d( 4)*le/2 + N20'*d(5) + N21'*d( 6)*le/2;
w(:,3) = N10'*d(5) + N11'*d( 6)*le/2 + N20'*d(7) + N21'*d( 8)*le/2;
w(:,4) = N10'*d(7) + N11'*d( 8)*le/2 + N20'*d(9) + N21'*d(10)*le/2;
w(:,5) = N10'*d(9) + N11'*d(10)*le/2 + N20'*d(11)+ N21'*d(12)*le/2;

% deflection overall beam, evaluated at 250 points
wg =-w(:);
dx = length/250;
x  = dx:dx:length;

plot(x,wg,'k-','LineWidth',4); 
fontsize = 18; set(gca,'FontSize',fontsize);
xlabel('x'); ylabel('deflection u_y (x)');
% print -depsc curve.eps
