function [ w ] = solveBeam( ~ )
%solveBeam.M gives deflection for the given boundary and loading conditions
% For ME256 3rd Homework
disp('This function is dependent on input files : node.dat, elem.dat, dispbc.dat, forces.dat')

load node.dat
load elem.dat
load forces.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)/12;
%
% 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
[U,R] = fembeam(A,L,E,nx,ny,ncon,NELEM,NNODE,F,dispID,Inertia);
w=U(2:3:end);
% plot(w,nx,'Linewidth',2)
end

