%% 1d linear elastic bar element
%
% Numerical Analysis of Structures WS 2017
%
% (c) 2017 Marek Tyburec, marek.tyburec@fsv.cvut.cz

%% DEFINE STRUCTURE
% nodal locations in global coordinates [m]
x = [2 1 0];

% nodal (external) forces in global coordinates [N]
nodeNumberFx = [2 3];
Fx = [-4 10];

% prescribed nodal displacements in global coordinates [m]
nodeNumberUx = [1];
Ux = [4];

% elements (node numbers)
elements = [3 2; 2 1];

% cross sectional areas of bars [m2]
A = [1 2];

% young modulus of bars [Pa]
E = [1 0.5];

% auxiliary variables
numberOfNodes = length(x);
numberOfElements = size(elements,1);

% prealocation
% number of degrees of freedom per node=1 (horizontal displacement)
stiffnessMatrixGlobal = zeros(1*numberOfNodes, 1*numberOfNodes);
elementLengths = zeros(numberOfElements,1);

% assembly of the global stiffness matrix
% for all elements
for element=1:numberOfElements
  % pick the element 'element'
  
  % find starting and ending node of the element
  FirstNode = elements(element, 1);
  SecondNode = elements(element, 2);
  
  % code numbers of the node: [FirstNode SecondNode]
  codeNumbersElement = [FirstNode SecondNode];
                        
  % length of the element
  elementLengths(element) = x(SecondNode) - x(FirstNode);
    
  % element stiffness matrix in global coordinates
  stiffnessMatrixElement = E(element) * A(element) / elementLengths(element) * ...
                              [ 1 -1; ...
                               -1  1];
                            
  % assembly (add element stiffness matrix to the structural global stiffness matrix)
  stiffnessMatrixGlobal(codeNumbersElement, codeNumbersElement) = ...
      stiffnessMatrixGlobal(codeNumbersElement, codeNumbersElement) + stiffnessMatrixElement;
  
end

% global force vector
forceVector = zeros(numberOfNodes,1);
forceVector(nodeNumberFx) = Fx;

% prescribed displacement vector
displacementVector = zeros(numberOfNodes,1);
displacementVector(nodeNumberUx) = Ux;

% essential code numbers (prescribed displacement)
essentialCodeNumbers = nodeNumberUx;

% degrees of freedom (F-nodes)
degreesOfFreedoms = setdiff(1:numberOfNodes, essentialCodeNumbers);

% solution
displacementVector(degreesOfFreedoms) = stiffnessMatrixGlobal(degreesOfFreedoms, degreesOfFreedoms) \ ... 
                                          (forceVector(degreesOfFreedoms) - ...
                                          stiffnessMatrixGlobal(degreesOfFreedoms,essentialCodeNumbers)*...
                                          displacementVector(essentialCodeNumbers))

% print unknown displacements to screen
for dof=1:length(degreesOfFreedoms)
  fprintf('u_%d = %f\n', degreesOfFreedoms(dof), displacementVector(degreesOfFreedoms(dof)));
end

reactions = stiffnessMatrixGlobal(essentialCodeNumbers, essentialCodeNumbers) * ...
            displacementVector(essentialCodeNumbers) + ...
            stiffnessMatrixGlobal(essentialCodeNumbers,degreesOfFreedoms) * displacementVector(degreesOfFreedoms);
            
for rea = 1:length(essentialCodeNumbers)
  fprintf('r_%d = %f\n', essentialCodeNumbers(rea), reactions(rea));
end
            
%% POST-PROCESSING
strains = (displacementVector(elements(:,2)) - displacementVector(elements(:,1)))./elementLengths;
for element = 1:numberOfElements
  fprintf('strain_%d = %f\n', element, strains(element));
end

stresses = E .* strains;
for element = 1:numberOfElements
  fprintf('stress_%d = %f\n', element, stresses(element));
end
