clear 
clc

% Solution of Assignement#1

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

%% Material properties
E = 70000;

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

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

%% Elements properties
areas = [645.2, 645.2, 645.2];  

element_data = {};  
%2dTruss: '2dTruss', nodes, coordNode1, coordNode2, E A
elem1 = {'2dTruss', connectivity(1, :), coordinates(connectivity(1, 1), :), coordinates(connectivity(1, 2), :), E, areas(1)};
elem2 = {'2dTruss', connectivity(2, :), coordinates(connectivity(2, 1), :), coordinates(connectivity(2, 2), :), E, areas(2)};
elem3 = {'2dTruss', connectivity(3, :), coordinates(connectivity(3, 1), :), coordinates(connectivity(3, 2), :), E, areas(3)};
element_data{1} = elem1;  
element_data{2} = elem2;
element_data{3} = elem3;

%% 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 = [1,2,4];    
freeDof = setdiff(allDof, fixedDof);
f_global(6) = -2000e3;           

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

%% Print the results
Rx_1=f_global(1)/1e3
Ry_1=f_global(2)/1e3
Ry_2=f_global(4)/1e3

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

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


test=1;