%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Finite Element MATLAB script for
% deflection and stress analysis of beams
% Written by: G.K. Ananthasuresh
% for MEAM 310 class in Spring 1998.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 clf					% Clear graphics window
 clear all				% Clear all variables
 clf					% Clear command window
 hold off				% No hold on the graphics window
%
% This script needs the following scripts to run
%	matcut.m	--> trims a matrix
%	veccut.m	--> trims a vector
%	femtruss.m	--> FEA for trusses
% It also needs the following input files	
%	node.dat	--> 	nodal data
%	elem.dat	-->	element data
%	forces.dat	-->	force data
%	disp.dat	-->	displacement boundary conditon data
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ INPUT from files 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Included global variables from file globals.m
%
% Read in nodal and element connectivity data from
% files: nodes.dat and elem.dat
load node.dat
load elem.dat
%
% Read in force data from the file forces.dat
load forces.dat
% Read in displacement boundary condition data from the file disp.dat
load dispbc.dat
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PRE-PROCESSING 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Identify the number of nodes, X and Y Coordinates of the nodes
NNODE = size(node,1);
nx = node(:,2); 
ny = node(:,3);
%
% Identify the number of elements and form an element connectivity array,
% the cross-section and Young's modulus arrays.
NELEM = size(elem,1);
ncon = elem(:,[2 3]);
beam_width = elem(:,4);
beam_depth = elem(:,5);
E = elem(:,6);
A = beam_width.*beam_depth;
Inertia = beam_width.*(beam_depth.^3);
%
% Arrange force information into a force vector, F
F = zeros(3*NNODE,1);			% Initialization
Nforce = size(forces,1);
for i = 1:Nforce,
   dof = (forces(i,2)-1)*3 + forces(i,3);
   F(dof) = forces(i,4);
end
%
% Displacement boundary conditions
Nfix = size(dispbc,1);
j = 0;
for i = 1:Nfix,
    j = j + 1;
	dispID(j) = (dispbc(i,2)-1)*3+dispbc(i,3);
        dispVal(j) = dispbc(i,4);
end
[dispID sortIndex] = sort(dispID);
dispVal = dispVal(sortIndex);
%
% Compute the lengths of the elements
for ie=1:NELEM,
  eye = ncon(ie,1);
  jay = ncon(ie,2);
  L(ie) = sqrt ( (nx(jay) - nx(eye))^2 + (ny(jay) - ny(eye))^2 );
end
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SOLUTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Call fembeam.m to solve for the following.
%	Deflections: U
%	Reaction forces at the constrained nodes: R
%	Internal forces in each truss member: Fint
%	Global stiffness matrix: Kglobal
%	Strain energy: SE
[U,R,Fint,Fint_local,Kglobal,SE] = fembeam(A,L,E,nx,ny,ncon,NELEM,NNODE,F, ...
   dispID,dispVal,Inertia); 
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% POST-PROCESSING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting
clf
% The "scale factor" below is used for scaling the displacements for the purpose
% of plotting. Change it as needed.
scale_factor = 10;
U = scale_factor * U;

for ip=1:NELEM,
  pt1 = ncon(ip,1); pt2 = ncon(ip,2);
  dx1 = U(3*(pt1-1)+1);
  dy1 = U(3*(pt1-1)+2);
  dx2 = U(3*(pt2-1)+1);
  dy2 = U(3*(pt2-1)+2);
% 
  plot([nx(pt1) nx(pt2)], [ny(pt1) ny(pt2)],':k');
  hold on
  plot([nx(pt1) nx(pt2)], [ny(pt1) ny(pt2)],'.k');
  plot([nx(pt1)+dx1 nx(pt2)+dx2], [ny(pt1)+dy1 ny(pt2)+dy2],'-r');
  plot([nx(pt1)+dx1 nx(pt2)+dx2], [ny(pt1)+dy1 ny(pt2)+dy2],'.r');
  hold on    
end
xlabel('X');
ylabel('Y');
axis('equal');

