% truss.m
%
% Matthias Heinkenschloss
% Sept 8, 2008

clear all
%
set(0, 'defaultaxesfontsize',16,'defaultaxeslinewidth',1,...
    'defaultlinelinewidth',1,'defaultpatchlinewidth',1)


% 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
figure(1); clf
% plot the bars for 
for i = 1:m
    x = [node(bar(i,1),1), node(bar(i,2),1)];
    y = [node(bar(i,1),2), node(bar(i,2),2)];
    plot(x,y,'b'); hold on
end
% plot the nodes
plot(node(:,1),node(:,2),'bo'); hold on



%
% 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
[ATKA] = stiff( bar, node, area, young, free);


% set right hand side
force     = zeros(2*n,1);
force(6)  = 100*9.80665;

% solve the linear system
x     = ATKA \ force(free);



% plot the deformed truss
% determine the locations of the node of the deformed truss
dnode = node;  % displaced nodes
mag   = 100;   % magnification facor for displacement
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) - mag*x(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) + mag*x(i);
    end
end

% plot the bars
 for i = 1:m
        x = [dnode(bar(i,1),1), dnode(bar(i,2),1)];
        y = [dnode(bar(i,1),2), dnode(bar(i,2),2)];
        plot(x,y,'r--'); hold on
end
% plot the nodes
plot(dnode(:,1),dnode(:,2),'ro'); hold on