%% FEM Example 2.1
% Sam Ciocco
clear, clc, close all;
%% Variable Inputs
numElements = 6; %Number of elements
numNodes = 5; %Number of nodes
DOF = 2; %Degrees of Freedom
A = 8; %Area in in^2
E = 1.9E6; %Elastic Modulus, in lb/in^2
forces = [0;0;0;0;0;0;0;-500;0;-500]; %Forces Vector, in order 1x;1y;2x;2y;...
constrained = [1,3]; %List of constrained nodes
%% Element Table
% For now, going to hardcode the table
% Will add node locations to calculate lengths and angles if I have time
elements = ones(numElements,7);
for i=1:numElements
elements(i,1) = i;
end
% Declare the first node for each element
elements(:,2) = [1;2;3;2;2;4];
% Declare the second node for each element
elements(:,3) = [2;3;4;4;5;5];
% Declare lengths for each element
elements(:,4) = [36;50.9;36;36;50.9;36];
% Declare angles for each element, in degrees
elements(:,5) = [0;135;0;90;45;0];
%Insert sine and cosine for each element, based on the angle
for j=1:numElements
elements(j,6) = cosd(elements(j,5));
elements(j,7) = sind(elements(j,5));
end
%Element table organized as follows:
% Element | Node 1 | Node 2 | Length (in) | Angle (deg) | Cosine | Sine
%% Set up each local matrix
%Apparently I figured out how to make an array of matrices
for i=1:numElements
keq(i) = A*E/elements(i,4);
Klocal(:,:,i) = keq(i).*[elements(i,6)^2,elements(i,6)*elements(i,7),-(elements(i,6)^2),-elements(i,6)*elements(i,7);
elements(i,6)*elements(i,7),elements(i,7)^2,-elements(i,6)*elements(i,7),-(elements(i,7)^2);
-(elements(i,6)^2),-elements(i,6)*elements(i,7),elements(i,6)^2,elements(i,6)*elements(i,7);
-elements(i,6)*elements(i,7),-(elements(i,7)^2),elements(i,6)*elements(i,7),elements(i,7)^2];
end
%% Set up Global Matrix
% Declare a matrix of correct size for each node to have its number of
% degrees of freedom (DOF)
Kglobal = zeros(numNodes*DOF,numNodes*DOF);
% Now painstakingly fill in each local matrix into the global matrix
% Ensuring each local matrix is placed correctly according to its nodes
for i=1:numElements
for j=0:1
for k=0:1
Kglobal(elements(i,2)*2-1+j,elements(i,2)*2-1+k) = Kglobal(elements(i,2)*2-1+j,elements(i,2)*2-1+k)+Klocal(j+1,k+1,i);
end
end
for j=0:1
for k=0:1
Kglobal(elements(i,2)*2-1+j,elements(i,3)*2-1+k) = Kglobal(elements(i,3)*2-1+j,elements(i,2)*2-1+k)+Klocal(j+1,k+3,i);
end
end
for j=0:1
for k=0:1
Kglobal(elements(i,3)*2-1+j,elements(i,3)*2-1+k) = Kglobal(elements(i,3)*2-1+j,elements(i,3)*2-1+k)+Klocal(j+3,k+3,i);
end
end
for j=0:1
for k=0:1
Kglobal(elements(i,3)*2-1+j,elements(i,2)*2-1+k) = Kglobal(elements(i,3)*2-1+j,elements(i,2)*2-1+k)+Klocal(j+3,k+1,i);
end
end
% Old code, not completely sure why it didn't work
% Kglobal(elements(i,2)*2-1,elements(i,2)*2-1) = Kglobal(elements(i,2)*2-1,elements(i,2)*2-1)+Klocal(1,1,i);
% Kglobal(elements(i,2)*2-1,elements(i,2)*2-1+1) = Kglobal(elements(i,2)*2-1,elements(i,2)*2-1+1)+Klocal(1,2,i);
% Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1) = Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1)+Klocal(2,1,i);
% Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1+1) = Kglobal(elements(i,2)*2-1+1,elements(i,2)*2-1+1)+Klocal(2,2,i);
end
%% Create list of unconstrained (aka moving) nodes
% This code assumes each node can either move or not; no partial
% constraints
nodes = 1:numNodes; %Create list of nodes
movingNum = 1;
moving = zeros(numNodes-length(constrained));
for i=1:numNodes
if((nodes(i) ~= constrained(1)) && (nodes(i) ~= constrained(2)))
moving(movingNum) = i;
movingNum = movingNum+1;
end
end
%% Trim Global Matrix in Accordance with Constrained Nodes
% Again assuming any constrained nodes are fully constrained
Kglobal2 = zeros((numNodes-length(constrained))*2,(numNodes-length(constrained))*2);
for i=1:length(moving)
for j=1:length(moving)
for k=0:1
for m=0:1
Kglobal2(i*2-1+k,j*2-1+m) = Kglobal(moving(i)*2-1+k,moving(j)*2-1+m);
end
end
end
end
% for l=1:length(moving)
% for i=:1
% for j=0:1
% Kglobal2(l*2-1+i,l*2-1+j) = Kglobal(moving(l)*2-1+i,moving(l)*2-1+j);
%
% end
% end
% end
%% Trim Forces in Accordance with Constrained Nodes
forces2 = zeros(length(moving)*2,1);
for l=1:length(moving)
forces2(l*2-1) = forces(moving(l)*2-1);
forces2(l*2) = forces(moving(l)*2);
end
%% Solve Matrix for movements
movements = Kglobal2\forces2;
%movements is the vector of the first node x, the first node y, the second
%node x, the second node y,...
numMoving = 1;
for i=1:2:length(movements)
fprintf('\nThe movement of node %.0f in the x-direction is %.3f inches\n',moving(numMoving),movements(i));
fprintf('The movement of node %.0f in the y-direction is %.3f inches\n',moving(numMoving),movements(i+1));
numMoving = numMoving+1;
end