%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Finite Element MATLAB script for
% deflection and stress analysis of beams
% Written by: G.K. Ananthasuresh
% for MEAM 310 class in Spring 1998 at the University of Pennsylvania
% Modified in August, 2014, for ME 237/NE211 at the Indian Institute of Science
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf					% Clear graphics window
clear all				% Clear all variables
clf					% Clear command window
hold off				% No hold on the graphics window
rho=2300;
%
% 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)/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
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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);
%

aa=forces(1,2); % node numner for calculating efective stiffness 
 
Uaa=U(3*(aa-1)+2); % displacement of node aa in y direction 
keff=forces(1,4)/Uaa
mm=0;
 
for ie=1:NELEM,
    if(elem(ie,9)==1)
       mm=mm+ U(3*(ie-1)+2)^2*L(ie)*A(ie);
    end
end
meff=rho*mm/Uaa^2



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 = 1;
Us = scale_factor * U;
for ip=1:NELEM,
    pt1 = ncon(ip,1); pt2 = ncon(ip,2);
    %
    uelem = [Us(3*(pt1-1)+1); Us(3*(pt1-1)+2); Us(3*pt1); Us(3*(pt2-1)+1); Us(3*(pt2-1)+2); Us(3*pt2)];
    plotbeamelem(nx(pt1),ny(pt1),nx(pt2),ny(pt2),uelem,20);
    hold on
end
xlabel('X');
ylabel('Y');
axis('equal');

