%{
  Deformation of a structure made of springs, all oriented and loaded
  in one direction. The structure is "solved" in terms of nodal
  DISPLACEMENTS w.r.t. the unloaded configuration.

  clear all; close all; clc;

  % The structure
  map = [1 2;2 3;2 3;3 4];
  kappa = [1000,1000,3000,1500];
  % The loads
  G = [0;0;5000;-1000];
  % The constraints [node ID, node displacement]
  C(1,1:2) = [1,0];
  % Solve
  [U,G,du]=Spring_Structure_1D(map,kappa,G,C);



%}
function [U,G,du]=Spring_Structure_1D(map,kappa,G,C)

%% Nr of elements and nodes
Nel = size(map,1);
Nn  = max(map(:));

%% Initialize solution vector and global stiffness matrix
U = nan(Nn,1);
K = sparse(Nn,Nn);

%% Assign boundary condition to constrained displacements
U(C(:,1)) = C(:,2);

%% Assembly
ONE = [1 -1;-1 1];
for e=1:Nel
    KE = ONE.*kappa(e);
    K(map(e,:),map(e,:)) = K(map(e,:),map(e,:)) + KE;
end

%% Enforce boundary conditions
IsNotConstrained = setdiff(1:Nn,C(:,1));

%% Solution
U(IsNotConstrained) = K(IsNotConstrained,IsNotConstrained)\G(IsNotConstrained);

%% Nodal reactions
G(C(:,1)) = K(C(:,1),:)*U;

%% Elongation of each spring
du = nan(Nel,1);
for e=1:Nel
    du(e) = U(map(e,2))-U(map(e,1));
end

end