%----------------------------------------------------------------------------
%                                 ENGR690
%                            
%
%                           Dr. Chris L. Mullen
%                                10-12-99
%
% Solve Generalized 1-D Model BVP in Elastodynamics:
%
%	a*u,xx + b*u,x + c*u = d*u,tt
%
%   id(ndof,numnp)  = boundary condition codes
%   f(ndof,numnp)   = prescribed nodal forces/displacements
%
%   lm(nee)         = localized element dofs                                         
%   s(nee,nee)      = element stiffness matrix 
%   m(nee,nee)      = element mass matrix
%   p(nee)          = element force vector
%
%   ak(neq,neq)      = global stiffness matrix 
%   am(neq,neq)      = global stiffness matrix
%   b(neq)          = global force vector                                                 
%           
%----------------------------------------------------------------------------            

%---------------------------------------------------------------------------
%  INPUT PHASE
%---------------------------------------------------------------------------
%-----------------------------------------
%  Problem Control parameters
%-----------------------------------------
numnp=5;              % total number of nodes in system
nsd=1;                % number of spatial dimensions
ndof=1;               % number of dofs per node
nel=4;                % number of elements
ned=1;                % number of element dofs per node
nen=2;                % number of element nodes
%-----------------------------------------
%  Nodal coordinate values
%-----------------------------------------
H=1
x(1)=0.0; x(2)=0.25; x(3)=0.5; x(4)=0.75; x(5)=1.0; 
x=x*H
%-----------------------------------------
%  Nodal connectivity for each element
%-----------------------------------------
ien(1,1)=1; ien(2,1)=2;
ien(1,2)=2; ien(2,2)=3;
ien(1,3)=3; ien(2,3)=4;
ien(1,4)=4; ien(2,4)=5;  
%-----------------------------------------
%  Wave velocity for case of shear beam
%-----------------------------------------
G=1;
ro=1;
A= 1;
k= 1;
Av=k*A
c2= k*G/ro
%-----------------------------------------
%  Coefficients of the ODE
%    u,xx= (1/c2)*u,tt
%-----------------------------------------
c=sqrt(c2)
aa=1; bb=0; cc=0; dd=1/c2
%-----------------------------------------
%  Prescribed nodal forces/boundary conditions
%-----------------------------------------
id=zeros(ndof,numnp);  % clear ID array
f=zeros(ndof,numnp);   % clear F array
%  (No nodal forces are prescribed for this problem)
id(1,1)=1;             % restrain dof 1 (base)
%-----------------------------------------
%  Nodal values of distributed force (pressure)
%-----------------------------------------
ff= zeros(ndof,numnp);
%---------------------------------------------------------------------------
%  INITIALIZATION PHASE
%---------------------------------------------------------------------------
gdof    = ndof*numnp;             % number of global dof
[id,neq]= bc(id,ndof,numnp);      % number of global equations
ak      = zeros(neq,neq);         % initialize global stiffness matrix
am      = zeros(neq,neq);         % initialize global stiffness matrix
[b]     = loads(f,id,ndof,numnp); % transfer nodal loads to global force 
%-----------------------------------------
%  Form     element arrays
%  Assemble global  arrays
%-----------------------------------------
nee=nen*ned;                        % number of element equations
for iel=1:nel                       % loop over all elements 
  xl   = lcoord(x,ien,nsd,nen,iel);    % localize nodal coordinates
  he   = xl(1,2)-xl(1,1);              % compute element length
  ffl  = lpress(ff,ien,ned,nen,iel);   % localize nodal pressure values
  [s,m,p]= el1d2(aa,bb,cc,dd,he,ffl'); % form element stiffness, mass,
				       %     and force
  s
  m
  lm   = local(id,ien,ned,nen,iel);    % localize ID array
  dl   = ldisbc(f,ien,lm,ned,nen,iel); % localize prescribed displacement bc
  al  =  zeros(nee)		       % localize prescribed acceleration bc
  p    = disbc1(p,s,m,dl,nee);         % add -"k"*g and -"m"*g,tt to force
  p
  [ak,am,b]= addstf1(ak,am,b,s,m,p,lm,nee); % assemble global stiffness, 
				       %      mass, and force
end
%---------------------------------------------------------------------------
%  SOLUTION PHASE
%---------------------------------------------------------------------------
%-----------------------------------------
%  FEM solution
%-----------------------------------------
ak
am
b'
[phi,w2]= eig(ak,am);             % solve eigenvalue problem (K-w2M)*phi=0
w=sqrt(w2)
freq=w/(2*pi);
w
freq=diag(freq)
fnfe(1)=freq(4)
fnfe(2)=freq(3)
fnfe(3)=freq(2)
fnfe(4)=freq(1)
phi
%f= disptr(d,f,id,ndof,numnp);   % transfer solution to F array
%-----------------------------------------
%  Exact solution
%-----------------------------------------
num=1:1:nel;
wn=(2*num-1)*(pi/2)*c/H
fn=(2*num-1)*(1/4)*c/H
diff=(fnfe-fn) 
%---------------------------------------------------------------------------
%  OUTPUT PHASE
%---------------------------------------------------------------------------
femvstheory=[num'  fnfe'    fn'   diff']
%---------------------------------------------------------------