%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%
% a beam model for testing eigenvalue analysis and 
% Programmed by K C Park, February 24, 2004
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%
%   Input data
%
% default data!!
L = 1.00;   % beam span
b = 0.05;   % beam width
h = 0.00;   % beam thickness
E = 69.e+9; % young's modulus
Iz = b*h^3/12; % second area moment  
EI = E*Iz;  % bending rigidity
V=b*h*L; rho = 2700;  m = rho*V;
node_g = 23; % default number of elements

disp('  ');
disp('Welcome to vibration analysis of a beam by the finite element method');
disp('Please provide the beam model data');
disp(' ');
L = input('beam span length (meter)  =  ');
b = input('beam width (meter)  =        ');
h = input('beam thickness (meter)  =    ');
E = input('Youngs modulus (Pascal)  =   ');
rho=input('density (kg/meter^3)  =      ');
nelem =input('number of elements to model the beam = ');
disp(' ')
disp('Please select the bounday conditions');
disp('free-free beam:   1');
disp('fixed-free beam:  2');
disp('fixed-fixed beam: 3');
disp('fixed-simply supported beam: 4');
disp('simple-simple beam: 5');
bccase =input('boundary condition forthis analysis  = ');
disp(' ');
n_mode =input('specify number of modes to becomputed  =  ');
%
% begin computations
%
node_g = nelem+1; % total number of global nodes
node_e = 2;   % number of nodes in each element
dof_node  = 2;   % degrees of freedom per node
dof_elem = dof_node*node_e;  % total elemental dofs
dof_g = dof_node*node_g;  % total global dofs
Iz = b*h^3/12; % second area moment  
EI = E*Iz;  % bending rigidity

%
%  element stiffness matrix
ell = L/nelem;  % elemental beam length
kr = EI/ell^3;
ke=kr*[
    12   6*ell         -12    6*ell;
   6*ell 4*ell^2       -6*ell 2*ell^2; 
   -12  -6*ell          12   -6*ell;
   6*ell 2*ell^2       -6*ell 4*ell^2];

mass=rho*b*h*ell;   
% lumped element mass matrix
me = sparse(zeros(dof_elem,dof_elem));
me(1,1)  = 0.5*mass;
me(3,3)  = 0.5*mass;
me(2,2)  = 0.5*mass*ell^2/12; me(4,4)  = 0.5*mass*ell^2/12;
%consistent mass case
me = mass*[
         13/35         11*ell/210        9/70          -13*ell/420;
         11*ell/210    ell^2/105        13*ell/420    -ell^2/140;
         9/70          13*ell/420        13/35         -11*ell/210;
         -13*ell/420   -ell^2/140        -11*ell/210   ell^2/105];	  
%
% generate connectivity table
%
conn=sparse(zeros(nelem, node_e));  % connectivity table 
%
IC=0;
for j = 1: node_g; 
IC = IC+1;
conn(IC, 1) = j;
conn(IC, 2) = j+1;
end;
% 
% generate assembly operator 
% 
Ls = sparse(zeros(nelem*dof_elem, dof_g));
i2 = speye(2);
for j = 1:nelem,
  for k = 1:node_e,
    r = (j-1)*dof_elem + (k-1)*dof_node; r = (r+1:r+dof_node);
    c = (conn(j,k)-1)*dof_node; c = (c+1:c+dof_node);
    Ls(r,c) = i2;
  end
end
%
% generate  the nodal coordinates
%
coord= zeros(node_g,1);
nodes = 0;
for i = 1: node_g  % x-direction sweep
  x =(i-1)*ell;
nodes = nodes+1;
coord(nodes, 1) = x;
end;

% obtain the unassembled stiffness matrices

K_ell = sparse(zeros(nelem*dof_elem, nelem*dof_elem));
M_ell = sparse(zeros(nelem*dof_elem, nelem*dof_elem));

for ne =1:nelem,
	  ig = (ne-1)*dof_elem +1: ne*dof_elem;
	  jg = (ne-1)*dof_elem +1: ne*dof_elem;
      K_ell(ig, jg) = K_ell(ig, jg) + ke;
	  M_ell(ig, jg) = M_ell(ig, jg) + me;
end;

% assemble global mass and stiffness matrices

Kg = Ls' * K_ell * Ls;
Mg = Ls' * M_ell * Ls;
%
%  apply boundary conditions - simply supported at both ends
%
ng = size(Kg,1);
if     (bccase==1)  % free free case
	wmodes_dof =1:2:(ng-1);
	BeamModel = 'Free - Free Beam';
	% do nothing
elseif (bccase==2)	% fixed - free case
	Kg([1,2],:) = []; Kg(:,[1,2]) = []; 
	Mg([1,2],:) = []; Mg(:,[1,2]) = []; 
	nn=size(Kg,1);
	wmodes_dof =1:2:(nn-1);
	BeamModel = 'Fixed - Free Beam';
elseif (bccase==3)	% fixed - fixed case
	Kg([1,2,(ng-1),ng],:) = []; Kg(:,[1,2,(ng-1),ng]) = []; 
	Mg([1,2,(ng-1),ng],:) = []; Mg(:,[1,2,(ng-1),ng]) = []; 
	nn=size(Kg,1);
	wmodes_dof =1:2:(nn-1);
	BeamModel = 'Fixed - Fixed Beam';
elseif (bccase==4)	% fixed - simple support case
	Kg([1,2,(ng-1)],:) = []; Kg(:,[1,2,(ng-1)]) = []; 
	Mg([1,2,(ng-1)],:) = []; Mg(:,[1,2,(ng-1)]) = []; 
	nn=size(Kg,1);
	wmodes_dof =1:2:(nn-1);
	BeamModel = 'Fixed - Simple Support Beam';
elseif (bccase==5)	% simple - simple support case
	Kg([1,(ng-1)],:) = []; Kg(:,[1,(ng-1)]) = []; 
	Mg([1,(ng-1)],:) = []; Mg(:,[1,(ng-1)]) = []; 
	nn=size(Kg,1);
	wmodes_dof =2:2:(nn-1);
	BeamModel = 'Simple - Simple Support Beam';
end;

%compute eigenvalues and eigenvectors

Minv = inv(Mg); MinvKg = Minv*Kg;

% obtain eigenvalues and eigenvectors . (non-mass normalized!)
if (bccase==1)
	% special case for free free beam
	[phi_total, lambda_total]= eig(full(Kg), full(Mg));
	 ndof = size(MinvKg,1);
	lambdaD  = diag(lambda_total);
	mode_extract = (ndof-n_mode+1):ndof;
	lambda=lambdaD(mode_extract);
	phi=phi_total(:,mode_extract);
	% convert to Hertz
	omega = sqrt(lambda);
	freq = omega/(2*pi);
	betaL = ((rho*b*h/EI)^(0.25))*sqrt(omega)*L; 
	%compute the equivalent beta*L

	%plot eigenvector, eigenvalues
	for nplot =1:n_mode,
	mode = n_mode+1-nplot;
	figure(nplot); 
	plot(phi(wmodes_dof,mode));
	xlabel('Beam  span');
	ylabel('Mode shape ');
	legend([ num2str(nplot), '-th mode ', ...
             ' beta*L = ', num2str(betaL(mode)), ...
             ' Freq(Hertz) = ', num2str(freq(mode))]);
	title([BeamModel]);
	end;

else	
	[phi, lambda, flag]= eigs(full(Kg), full(Mg),n_mode,0);
	lambda = diag(lambda);
	% convert to Hertz
	omega = sqrt(lambda);
	freq = omega/(2*pi);
	betaL = ((rho*b*h/EI)^(0.25))*sqrt(omega)*L; 
	%compute the equivalent beta*L

	%plot eigenvector, eigenvalues
	for nplot =1:n_mode,
	figure(nplot); 
	plot(phi(wmodes_dof,nplot));
	xlabel('Beam  span');
	ylabel('Mode shape ');
	legend([ num2str(nplot), '-th mode ', ...
             ' beta*L = ', num2str(betaL(nplot)), ...
             ' Freq(Hertz) = ', num2str(freq(nplot))]);
	title([BeamModel]);
	end;
end; %end if loop