bridge.m

LU decomposition routines in matlab

% bridge
% set up and solve the equations for the bridge-truss 
%
% Matthias Heinkenschloss
% Jan 26, 2001

% set number of bar (m) and number of node (n)
m = 13;
n = 8;


% 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   3;
    1   2;
    2   3;
    3   4;
    5   3;
    5   2;
    5   4;
    4   7;
    5   7;
    5   6;
    6   7;
    8   7;
    6   8];
   
% 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]

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


% plot the undeformed and the deformed truss
truss_plot(bar, node, 1)
%print -depsc bridge.eps

% area is an m array. The i-th element contains the cross sectional
% area of bar i. We assume that all bars have rectangular cross
% section of size 0.01m times 0.05m, i.e., area = 0.0005 [m^2].
area = 0.0005*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 15 16];

% 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;

% right hand side
b     = zeros(2*n,1);
b(4)  = 1e6;
b(10) = 1e6;
b(12) = 1e6;


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)
    k = free(i);
    if( mod(k,2) == 0 ) 
        % k is even, i.e., u(i) is vertical displacement of node k/2
        dnode(k/2,2) = node(k/2,2) - u(i);
    else
        % k is odd, i.e., u(i) is horizontal displacement of node (k+1)/2
        dnode((k+1)/2,1) = node((k+1)/2,1) + u(i);
    end
end

% plot the undeformed and the deformed truss
truss_plot(bar, node, 1, dnode)
%print -depsc bridge_def.eps