% Introduction of Finite Elements in Engineering Mechanics 
% Dr. X. Frank Xu, Stevens Institute of Technology
% fe_bar1.m 


%----------------------
%Mesh Information
%----------------------
NumNP = 3; %number of nodal points
NumEL = 2; % number of elements
NP=[1 2;2 3]; %Nodes in elements from 1 to NumEL
XORD(1:NumNP,1)=[0,1,2]'; % X-coordinate of nodal points

%----------------------
%Specify Boundary Conditions
%----------------------
U=zeros(NumNP,1);
F=zeros(NumNP,1);
NPBC=zeros(NumNP,1);
NPBC(1)=1; %Displacement of Node 1 is constrained
U(1)=0; % Specified displacement for Node 1 
F(3)=1000; % Specified load for Node 3 (N)

%----------------------
%Specify Material Properties
%----------------------
E(1:NumEL,1)=1e11*[1,1]'; % Young's modulus (Pa)
A(1:NumEL,1)=1e-4*[2,1]'; % Cross section area(m2) 
L(1:NumEL,1)=1e-1*[1,1]'; % Length (m)

%----------------------
% Form Stiffness Matrix
%----------------------
K=zeros(NumNP, NumNP);
for I=1:NumEL
    Ke=zeros(2,2);
    k=E(I)*A(I)/L(I); % Bar stiffness
    Ke=[k -k;-k k];  %Element stiffness matrix
    PosK=[NP(I,1);NP(I,2)]; % Position of Ke in global stiffness matrix K 
    K(PosK,PosK)=K(PosK,PosK)+Ke; % Global matrix 
end
K0=K;
%--------------------------------------
% IMPLEMENTATION OF BOUNDARY CONDITIONS
%--------------------------------------
for I=1:NumNP
    if NPBC(I)==1 
       K(I,I)=K(I,I)*10000000; % Big number approach
       F(I)=K(I,I)*U(I);
    end
end

%----------------------
% Solve displacement
%----------------------
U = K\F;

%-----------------------------------------
% Post-processing to obtain stress and reaction forces
%-----------------------------------------
for I=1:NumEL
    def(I)=U(NP(I,2))-U(NP(I,1)); % Deformation of element I
    stress(I)=E(I)*def(I)/L(I); % Stress of element I
    force(I)=stress(I)*A(I); % Force of element I
end
F=K0*U; % Boundary forces