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

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

% nodal (external) forces in global coordinates [N]
nodeNumberFx = [2];
nodeNumberFy = [];
Fx = [10000];
Fy = [];

% prescribed nodal displacements in global coordinates [m]
nodeNumberUx = [1 3];
Ux = [0 0];
nodeNumberUy = [1 3];
Uy = [0 0];

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

% cross sectional areas of bars [m2]
A = [5*10^-4 5*10^-4];

% young modulus of bars [Pa]
E = [200*10^9 200*10^9];

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

% prealocation
% number of degrees of freedom per node=2 (two displacements)
stiffnessMatrixGlobal = zeros(2*numberOfNodes, 2*numberOfNodes);
cosinusVector = zeros(numberOfElements,1);
sinusVector = zeros(numberOfElements,1);
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_x, FirstNode_y, SecondNode_x, SecondNode_y]
  codeNumbersElement = [FirstNode*2-1 FirstNode*2 SecondNode*2-1 SecondNode*2];
                        
  % differences and length of the element
  deltaX = x(SecondNode) - x(FirstNode);
  deltaY = y(SecondNode) - y(FirstNode);
  elementLengths(element) = sqrt(deltaX^2 + deltaY^2);
  
  % cosinus and sinus, store it for post-processing
  cosinusVector(element) = deltaX/elementLengths(element);
  sinusVector(element) = deltaY/elementLengths(element);
  c = cosinusVector(element);
  s = sinusVector(element);
  
  % element stiffness matrix in global coordinates
  stiffnessMatrixElement = E(element) * A(element) / elementLengths(element) * ...
                              [ c^2  c*s -c^2 -c*s; ...
                                c*s  s^2 -c*s -s^2;
                               -c^2 -c*s  c^2  c*s;
                               -c*s -s^2  c*s  s^2];
                            
  % 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*2,1);
forceVector(2*nodeNumberFx-1) = Fx;
forceVector(2*nodeNumberFy) = Fy;

% prescribed displacement vector
displacementVector = zeros(numberOfNodes*2,1);
displacementVector(2*nodeNumberUx-1) = Ux;
displacementVector(2*nodeNumberUy) = Uy;

% essential code numbers (prescribed displacement)
essentialCodeNumbers = sort([nodeNumberUx*2-1 nodeNumberUy*2]);

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

% solution
displacementVector(degreesOfFreedoms) = stiffnessMatrixGlobal(degreesOfFreedoms, degreesOfFreedoms) \ ... 
                                          (forceVector(degreesOfFreedoms) - ...
                                          stiffnessMatrixGlobal(degreesOfFreedoms,essentialCodeNumbers)*...
                                          displacementVector(essentialCodeNumbers))
                                          
% plotting (uncommenting the following lines will plot the elements
%hold on;
%plot(x(elements)', y(elements)', 'Color', 'k', 'LineWidth', 2);
%scalingFactor = 100;
%displacedX = x + scalingFactor*displacementVector(1:2:end)';
%displacedY = y + scalingFactor*displacementVector(2:2:end)';
%plot(displacedX(elements)', displacedY(elements)', 'Color', 'r', 'LineWidth', 1, 'LineStyle','--');
%axis equal off;
%hold off;
