% truss.m
% set up and solve the equations for the truss in figure 1.3
% of section 1.9
%
% Matthias Heinkenschloss
% Jan 26, 2001


% set number of bar (m) and number of node (n)
m = 5;
n = 4;

% bar is an m times 2 array. The i-th row contains
% information for bar i in the following form:
%  [ lower-left-node   upper-right-node ]

bar ...
= [ 1   2;
    3   2;
    4   3;
    3   1;
    4   2];
   
% node is an n times 2 array. The i-th row contains
% information for node i in the following form:
%  [x-coordinate-of-node   y-coordinate-of-node]
% Coordinates are given in meters so that our units
% are compatible.

node ...
= [ 0   1;
    2   1;
    2   0;
    0   0];

% plot the undeformed truss
truss_plot(bar, node, 1)


% area is an m array. The i-th element contains the cross sectional
% area of bar i. We assume that all bars have circular
% cross section with diameter 0.5cm, i.e., area = 0.0025^2*pi [m^2].
area = (0.0025^2*pi)*ones(m,1);

% young is an m array. The i-th element contains the Young's
% modulus for bar i. We assume that all bars are made of the same 
% material with Young's modulus is 195 GPa 
young = 195.e9*ones(m,1);     


% set indices of fixed displacements
fixed = [1 2 7 8 ];

% determine indices of free displacements
free = [];
for i = 1:2*n
    if (~any(fixed==i))
       free = [free; i];
    end
end  

% determine the stiffness matrix
[K] = stiff( bar, node, area, young, free);

% factor the stiffness matrix
[K, ipivt, iflag] = lu_pp( K );


disp('Hit return to apply a load.......'); pause;
% set right hand side
b     = zeros(2*n,1);
b(6)  = 100*9.80665;

% solve the linear system
u                 = b(free);
[u, iflag]        = lu_pp_sl( K, u, ipivt );
if( iflag ~= 0 )
    error( ['    lu_pp_sl returned with iflag = ', int2str(iflag)])
end

% determine the locations of the node of the deformed truss
dnode = node;
for i = 1:size(free(:),1)
    j = free(i);
    if( mod(j,2) == 0 ) 
        % j is even, i.e., u(i) is vertical displacement of node j/2
        dnode(j/2,2) = node(j/2,2) - u(i);
    else
        % j is odd, i.e., u(i) is horizontal displacement of node (j+1)/2
        dnode((j+1)/2,1) = node((j+1)/2,1) + u(i);
    end
end

% plot the undeformed and the deformed truss
truss_plot(bar, node, 1, dnode)