function [d,F,K,k11,f1] = SolveLinearSystem(Node,Element,Load,p)
% [d,F,K,k,f] = SolveLinearSystem(Node,Element,Load,p)

% Walter Ponge-Ferreira
% Cotia, 17.06.2018

if nargin<4,
  p = 0;
end

n = GetNodeNumber(Node);
m = GetElementNumber(Element);
r = GetLoadNumber(Load);

if p == 1,
  % Print Nodes
  fprintf('\nNodes:\n\n');
  fprintf('No. \t X \t Y \t Cy \t Cz \n');
  for i = 1:n,
    fprintf('%d \t %g \t %g \t %d \t %d \n',i,Node(i,1),Node(i,2),Node(i,3),Node(i,4));
  end
  fprintf('\n');
  % Print Elements
  fprintf('\nElements: \n\n');
  fprintf('No. \t I1 \t I2 \t Type \t Par1 \t Par2 \t Par3 \n');
  for i = 1:m,
    fprintf('%d \t %d \t %d \t %s \t',i,Element{1,i}, Element{2,i}, Element{3,i});
    fprintf('%g \t %g \t %g \n',Element{4,i},Element{5,i},Element{6,i});
  end
  fprintf('\n');
  % Print Loads
  fprintf('\nLoads: \n\n');
  fprintf('No. \t Node \t Dir \t Value \n');
  for i = 1:r,
    fprintf('%d \t %d \t %d \t %g \n',i,Load(i,1),Load(i,2),Load(i,3));
  end
  fprintf('\n');
end

% Assemble Global Stiffness Matrix
K = StiffnessMatrixAssemble(Node,Element);

% Boundary Conditions
BC = reshape(Node(:,3:4)',1,2*n);

% Split the Stiffness Matrix K = [ k11 k12; k21 k22]
I = find(BC == 0);
J = find(BC != 0);

k11 = K(I,I);
k12 = K(I,J);
k21 = K(J,I);
k22 = K(J,J);

% Read Load Vector
F = zeros(2*n,1);
d = zeros(2*n,1);

for i = 1:r,
  if (Load(i,4) == 0),
    F(2*Load(i,1)+Load(i,2)-3)=Load(i,3);
  elseif (Load(i,4) == 1),
    d(2*Load(i,1)+Load(i,2)-3)=Load(i,3);
  else
    error('Incorrect load type');
  endif
end 

% Split Load Vector
f1 = F(I);

% Split Displacement Vector
d2 = d(J);

% Solve Linear Sistem
d1 = k11\(f1-k12*d2);
f2 = k21*d1 + k22*d2;

d(I) = d1;

% Applied + Reaction Forces
F = K*d;

endfunction