FEM 2D Truss Problem

% MATLAB Code for 2D Truss to find displacement and reactions using FEM
% =====================================================================
clear;
clc;

% Example: 2
% ===========
ND = [1 2 3 4 5]; % Node Number
NE = [1 2 3 4 5 6 7]; % Element Number
NC = [1 2;2 3;3 4;4 5;1 4;2 5;2 4]; % Node Connectivity
EL = [1 1 sqrt(2) 1 sqrt(2) sqrt(2) 1]; % Member Length (m)
alf = [0 0 3*pi/4 pi pi/4 3*pi/4 pi/2]; % Member Angle
E = 2.1E11; % Youngs Modulus (N/m2)
A = pi*(12.5/1000)^2/4; % Crossectional Ares (m2)
BC = [1 2 4 9 10]; % Boundary Condition
F_val = [0 200 -100 0 -100]; % Applied Load (kN)

n_D = length(ND);
n_E = length(NE);
dof = n_D*2;
dof_n = [(2.*ND-1)' (2.*ND)'];
KG = zeros(dof,dof);

for ii = 1:n_E
n1 = NC(ii,1);
n2 = NC(ii,2);
loc = [dof_n(n1,:) dof_n(n2,:)];
nj = length(loc);
c = cos(alf(ii));
s = sin(alf(ii));
T = [c s 0 0;-s c 0 0;0 0 c s;0 0 -s c];
ke = (E*A/EL(ii)).*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0];
Ke = T'*ke*T;
for jj = 1:nj
l1 = loc(jj);
for kk = 1:nj
l2 = loc(kk);
KG(l1,l2) = KG(l1,l2)+Ke(jj,kk);
end
end
end

K_temp = KG;
for jj = 1:length(BC)
K_temp(:,BC(jj)) = 0;
K_temp(BC(jj),:) = 0;
end
K_temp(~any(K_temp,2),:) = [];
K_temp(:,~any(K_temp,1)) = [];
KR = K_temp;

FR = F_val';
UR = KR\(FR.*1000);
UR = UR.*1000;
UG = [0 0 UR(1) 0 UR(2) UR(3) UR(4) UR(5) 0 0]';
Rc = (KG*(UG/1000))/1000;
disp('Displacement in mm')
disp(UR)
disp('Reaction in kN')
disp(Rc)