%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This m-file forms the finite element matrix equation of a bar and solves it.
% G.K. Ananthasuresh for ME 256 class in Jan.-May semester of 2011.
%
% Inputs:
%			Ae	--> A vector containing the cross-section areas of all the elements
%			Le	-->	A vector containing the lengths of all the elements
%			Ee	-->	A vector containing the Young's moduli of all the elements
%			nx	-->	A vector containing the X-coordinates of all the nodes
%			ny	-->	A vector containing the Y-coordinates of all the nodes
%		  ncon	-->	A matrix containing the nodal-connectivity of all the elements
%			F	-->	A vector containing the applied external forces
%		dispID	--> A vector containing the disp. boundary conditions
%
% Outputs:
%		u		-->	The displacements of all the dof
%		Rdisp	-->	The reaction forces at the dof which are specified to be fixed
%		P		--> A matrix containing theinternal reaction forces at the two
%					nodes of each element, and flag that tells you if the element
%					is in tension (1), compression (-1) or no stress (0).
%		King	--> The global stiffness matrix before the application of the
%					boundary conditions.
%		SE		--> The strain energy in the entire truss
%
function [u,Rdisp,P,Ksing,SE] = fembar(Ae, Le, Ee, nx, ny, ncon, NELEM, ...
    NNODE, F, dispID,dispVal)

K = zeros(NNODE,NNODE);		% Initialize the global stiffness matrix
k = zeros(2,2);				% Initialize local stiffness matrix

for ie=1:NELEM,
    eye = ncon(ie,1);
    jay = ncon(ie,2);
    % Form the local element matrix
    E = Ee(ie);
    A = Ae(ie);
    L = Le(ie);
    k(1,1)=E*A/L;
    k(1,2) = -k(1,1);
    k(2,1) = k(1,2);
    k(2,2) = k(1,1);
    % Form ID matrix to assemble klocal into the global stiffness matrix, K.
    id1 = eye
    id2 = jay
    K(id1,id1) = K(id1,id1) + k(1,1);
    K(id1,id2) = K(id1,id2) + k(1,2);
    K(id2,id1) = K(id2,id1) + k(2,1);
    K(id2,id2) = K(id2,id2) + k(2,2);
end

% Store unlaltered K as Ksing before applying the boundary conditions.
Ksing = K;

%det(K)
%inv(K);
%pause

% Imposing displacement boundary conditions
% ------------------------------------------
% dispID array contains the dof which are assigned specified values.
[sm,sn] = size(dispID);
Ndbc = sn;
for nd=1:Ndbc,
    for nr=1:NNODE-nd+1,
        F(nr) = F(nr) - K(nr,dispID(nd)-nd+1) * dispVal(nd);
    end
    K = matcut(K,dispID(nd)-nd+1);
    F = veccut(F,dispID(nd)-nd+1);
end

% To solve for unknown displacements.
U = inv(K)*F;
SE = .5*U'*K*U;

% Results
% ---------------
% "u" for all nodes (including those where values were specified)
u = zeros(NNODE,1);

for iu=1:Ndbc,
    u(dispID(iu)) = 12345.12345;
end

iuc = 0;
for iu=1:NNODE,
    if u(iu) == 12345.12345
        iuc = iuc+1;
    else
        u(iu) = U(iu-iuc);
    end
end

for iu=1:Ndbc,
    u(dispID(iu)) = dispVal(iu);
end


%----------------------------------------------
% Computation of reactions at constrained nodes
%----------------------------------------------
R = Ksing*u;
Rdisp = zeros(1,Ndbc);
for iu=1:Ndbc,
    Rdisp(iu) = R(dispID(iu));
end

%-------------------------------------------
% Computation of internal reaction forces
% and storing in P(NNODE,4)
%-------------------------------------------
%M = zeros(NNODE,1);

for ie=1:NELEM,
    
    eye = ncon(ie,1);
    jay = ncon(ie,2);
    
    for ie=1:NELEM,
        eye = ncon(ie,1);
        jay = ncon(ie,2);
        % Form the local element matrix
        E = Ee(ie);
        A = Ae(ie);
        L = Le(ie);
        k(1,1)=E*A/L;
        k(1,2) = -k(1,1);
        k(2,1) = k(1,2);
        k(2,2) = k(1,1);
        % Form ID matrix to assemble klocal into the global stiffness matrix, K.
        id1 = eye;
        id2 = jay;
        uelem = [u(id1) u(id2)];
        
        Rint = kelem*uelem';
        
        P(ie,1) = Rint(1);
        P(ie,2) = Rint(2);
    end  
end   
    