% Introduction of Finite Elements in Engineering Mechanics 
% Dr. X. Frank Xu, Stevens Institute of Technology
% femem_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 NUMNP
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 (unit: m)
F(3)=1000; % Specified load for Node 3 (unit: N)

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

%----------------------
% Form Stiffness Matrix
%----------------------
K=zeros(NumNP, NumNP);
for I=1:NumEL
    Ke=zeros(2,2);
    L(I)=XORD(NP(I,2))-XORD(NP(I,1)); % Length (unit: m)
    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

%--------------------------------------
% 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;