% Introduction of Finite Elements in Engineering Mechanics 
% Dr. X. Frank Xu, Stevens Institute of Technology
% simple beam theory
clear
%----------------------
%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,2,4]'; % X-coordinate of nodal points

%----------------------
%Specify Boundary Conditions
%----------------------
NumEQ = 2*NumNP;
U=zeros(NumEQ,1);
F=zeros(NumEQ,1);
NPBC=zeros(NumEQ,1);
NPBC(1)=1; %Y-Displacement of Node 1 is constrained
NPBC(2*3-1)=1; %Y-Displacement of Node 3 is constrained
U(1)=0; % Specified X-displacement for Node 1 
U(5)=0; % Specified Y-displacement for Node 3 

F(2*2-1)=-10e3; % Specified Y-force for Node 2 (unit:N)

%----------------------
%Specify Material Properties
%----------------------
E(1:NumEL,1)=210e9; % Young's modulus (unit:Pa)
II(1:NumEL,1)=.26e-6; % Moment inertia of cross-section(unit:m4) 

%----------------------
% Form Stiffness Matrix
%----------------------
K=zeros(NumEQ, NumEQ); 
for I=1:NumEL
    Ke=zeros(4,4);
    L=abs(XORD(NP(I,2))-XORD(NP(I,1)));
    Ke=II(I)*E(I)/L^3*[12 6*L -12 6*L;6*L 4*L^2 -6*L 2*L^2;-12 -6*L 12 -6*L;6*L 2*L^2 -6*L 4*L^2]; %Element stiffness matrix
    for J=1:2  
        PosK(J*2-1)=NP(I,J)*2-1;% Position of Ke in global stiffness matrix K 
        PosK(2*J)=NP(I,J)*2; % Position of Ke in global stiffness matrix K 
    end
    K(PosK,PosK)=K(PosK,PosK)+Ke; % Global matrix 
end
K0=K;
%--------------------------------------
% IMPLEMENTATION OF BOUNDARY CONDITIONS
%--------------------------------------
for I=1:NumEQ
    if NPBC(I)==1 
       K(I,I)=K(I,I)*1E9; % 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
%-----------------------------------------

FBC=K0*U; % Boundary forces  
       
beam_plot