clear 
clc

% Solution of Exercice#3

%% Dimension and degrees of freedom of the problem
nDim = 2;
nDof = 3; % nDof per node
nNodesPerElement = 2;

%% Material properties
E = 200000;

%% Node coordinates
ptA = [0, 0] * 1e3;
ptB = [2, 2] * 1e3;
ptC = [-4, 2] * 1e3;
coordinates = [ptA; ptB; ptC];
nNodes = size(coordinates, 1);
nDofTot=nNodes*nDof;

%% Connectivity
element1 = [1, 2]; 
element2 = [1, 3];
element3 = [1, 1];
element4 = [2, 2];
element5 = [3, 3];
connectivity = [element1; element2; element3; element4; element5];
nElements = size(connectivity, 1);

%% Elements properties
areas = [1000, 4000];
inertias = [3000, 12000];
kSprings = [50000.0, 1000000.0, 1000000.0];

element_data = {}; % Cell array to contain the properties of each element
% '2dTruss': {'2dTruss', nodes, coordNode1, coordNode2, E, A}
% '2dBeam': {'2dBeam', nodes, coordNode1, coordNode2, E, A, I}
% 'spring': {'spring', node, globalDof, k} where globalDof=1 for X, 2 for Y, 3 for rotation
elem1 = {'2dBeam', connectivity(1,:), coordinates(connectivity(1,1),:), coordinates(connectivity(1,2),:), E, areas(1), inertias(1)};
elem2 = {'2dBeam', connectivity(2,:), coordinates(connectivity(2,1),:), coordinates(connectivity(2,2),:), E, areas(2), inertias(2)};
elem3 = {'spring', connectivity(3,:), 2, kSprings(1)};
elem4 = {'spring', connectivity(4,:), 3, kSprings(2)};
elem5 = {'spring', connectivity(5,:), 3, kSprings(3)};

% Append elements to element_data
element_data{end+1} = elem1;
element_data{end+1} = elem2;
element_data{end+1} = elem3;
element_data{end+1} = elem4;
element_data{end+1} = elem5;


%% Equation number matrix
numEquations=computeEquationNumber(nDof,element_data);

%% Assemble global stiffness matrix
K_global=assembleGlobalStiffnessMatrix(nDofTot,numEquations,element_data);

%% Assemble load and displacement vectors
f_global = zeros(nDofTot, 1);
u_global = zeros(nDofTot, 1);

%% Apply boundary conditions
allDof = 1:(nDofTot); 
fixedDof = [4, 5, 7, 8];    
freeDof = setdiff(allDof, fixedDof);
u_global(2) = -5;           

% Solve the system to get forces and displacements
f_global(freeDof) = K_global(freeDof,freeDof) * u_global(freeDof);

%% Print the results
f_x=f_global(1)/1e3
f_y=f_global(2)/1e3
anglAlpha=atan(f_global(2)/f_global(1))*180/pi

%% Plot deformed structure
% plotDeformedStructure2d(coordinates,connectivity,u_global,1e2,'deformedStructure');

%% Plot internal forces
%plotInternalForces2dTruss(coordinates,connectivity, E, areas, u_global, 'axialLoad');


test=1;