% define coordinates of each nodes
Coord_nodes = [0,0; 10,0; 10,10; 0,10];
% define connection of each element
Connect_elements = [1,2; 1,3; 1,4; 4,3; 4,2; 2,3];
% define degree of Freedom of each nodes
Deg_free = [1,2; 6,5; 7,8; 3,4]; % [x-direction DOF, y-direction DOF]
% define number of reactions at supports
Nr = 3;
% define number of members
Nm = size(Connect_elements,1);
% define number of nodes
Nn = 4;
% define EA
EA = 1;
% creation of ID Array
ID = zeros(Nm,4);
for k = 1:Nm
i = Connect_elements(k,1);
j = Connect_elements(k,2);
ID(k,1:2) = Deg_free(i,1:2);
ID(k,3:4) = Deg_free(j,1:2);
end
% assemble Global Stiffness matrix
k_elem = nan(4,4,Nm);
for k = 1:Nm
i = Connect_elements(k,1);
j = Connect_elements(k,2);
dx = Coord_nodes(j,1) - Coord_nodes(i,1);
dy = Coord_nodes(j,2) - Coord_nodes(i,2);
L = norm([dx,dy]);
Lambda_x = dx/L;
Lambda_y = dy/L;
%
k11 = Lambda_x^2/L;
k12 = (Lambda_x*Lambda_y)/L;
k13 = -k11;
k14 = -k12; % not (-Lambda_x*Lambda_y)/L;
k21 = k12;
k22 = Lambda_y^2/L;
k23 = -k12;
k24 = -k22;
k31 = -k11;
k32 = k23;
k33 = k11;
k34 = k12;
k41 = k14;
k42 = k24;
k43 = k12;
k44 = k22;
k_elem(:,:,k) = EA * [...
k11, k21, k31, k41;...
k12, k22, k32, k42;...
k13, k23, k33, k43;...
k14, k24, k34, k44];
end
%
% Get total number of DOFs
Nd = max(Deg_free(:));
% Initialize global stiffness matrix
GSM = zeros(Nd,Nd);
% Assemble global stiffness matrix
for k = 1:Nm
for i = 1:4
for j = 1:4
GSM(ID(k,i), ID(k,j)) = GSM(ID(k,i), ID(k,j)) + k_elem(i,j,k);
end
end
end
Giving
display(GSM)
