clear 
clc

% Solution of Exercice#2

%% Dimension and degrees of freedom of the problem
nDim = 2;
nDof = 2; % 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];
connectivity = [element1; element2];
nElements = size(connectivity, 1);

%% Elements properties
areas = [1000, 4000];  

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)};
element_data{1} = elem1;  
element_data{2} = elem2;

%% 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 = [3, 4, 5, 6];    
freeDof = setdiff(allDof, fixedDof);
u_global(2) = -5;           

% Solve the system to get forces and displacements
f_global = K_global * u_global;

%% 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;