% Workshop Example 08  (WorkshopEx08.m)
%
% Build and Analyze a Determinate Truss
%
% This script evaluates a file that specifies a truss.
% It does this by defining the following four matrices
% 
% nodes - an nn-by-2 matrix with one row per node, where
%         each row specifies the (x,y) coordinates of the
%         corresponding node
%
% bars  - an nb-by-2 matrix with one row per bar, where
%         each row specifies the indices of the source and 
%         target nodes for the corresponding bar.
%
% fixed - an nn-by-2 matrix with one row per node, where
%         each row specifies whether there is a reactive
%         force in the horizontal and vertical directions.
%         [0 0]: If row i consists of [0 0] then node i has no
%                reactive forces.
%         [0 1]: If row i consists of [0 1] then node i has a
%                vertical reactive force and no horizontal reactive
%                force.  Node i can move horizontally but not
%                vertically.
%         [1 0]: If row i consists of [1 0] then node i has a
%                horizontal reactive force and no vertical reactive
%                force.  Node i can move vertically but not
%                horizontally.
%         [1 1]: If row i consists of [1 1] then node i has 
%                reactive forces in both the horizontal and vertical
%                directions.  Node i is immovable.
%
% loads - an nn-by-2 matrix with one row per node, where
%         each row specifies the horizontal and vertical loads
%         applied to that node.  
%         If row i consists of [h v] then node i has a horizontal
%         applied load of h and a vertical applied load of v.
%         Both h and v can be 0.
%
% The script then
% 1.  plots the specified truss
% 2.  builds the 2*nn-by-2*nn matrix that defines the sum of forces
%     for each node. (one equation for the horizontal components
%     and one equation for the vertical components)
% 3.  builds the vector of applied loads
% 4.  solves the system of equation for the internal and reactive
%     forces.

% Daniel D. Warner
% July 9, 2008

name = input('Enter the name of the file defining the truss: ','s');
eval(name)

% Plot the nodes
[nn,two] = size(nodes);
xn = nodes(:,1);
yn = nodes(:,2);
plot(xn,yn,'ok')

title(name)
% Set the view, include a small margin around the truss 
xmargin = 0.1*max(abs(xn));
ymargin = 0.1*max(abs(yn));
xmin = min(xn) - xmargin;
xmax = max(xn) + xmargin;
ymin = min(yn) - ymargin;
ymax = max(yn) + ymargin;
axis([xmin xmax ymin ymax]);
hold on

% Plot the bars
[nb,two] = size(bars);
for k=1:nb
    % Get the index for the first node (source node) of bar k
    sindx = bars(k,1);
    % Get the index for the second node (target node) of bar k
    tindx = bars(k,2);
    plot([xn(sindx) xn(tindx)],[yn(sindx) yn(tindx)])
end

% Label the nodes
for k=1:nn
    text(xn(k)+0.2*xmargin, yn(k)+0.2*ymargin, num2str(k))
end
hold off
pause

% Set up the system of 2*nn equations
A = zeros(2*nn,2*nn);
for k=1:nb
    % Determine the indices for the source and target nodes
    sindx = bars(k,1);
    tindx = bars(k,2);
    % Calculate the length of the bar
    leng = sqrt((xn(tindx)-xn(sindx))^2 + (yn(tindx)-yn(sindx))^2);
    % Calculate ck and sk
    ck = (xn(tindx)-xn(sindx))/leng;
    sk = (yn(tindx)-yn(sindx))/leng;
    % Add in the contributions to the horizontal and
    % vertical sums for the source node
    A(2*sindx-1,k) = ck;
    A(2*sindx,k) = sk;
    % Add in the contributions to the horizontal and
    % vertical sums for the target node
    A(2*tindx-1,k) = -ck;
    A(2*tindx,k) = -sk;    
end

% Add in the reactive forces for the fixed nodes
index = nn+1;
for k=1:nn
    if (fixed(k,1) == 1)
        A(2*k-1,index) = 1;
        index = index+1;
    end
    if (fixed(k,2) == 1)
        A(2*k,index) = 1;
        index = index+1;
    end
end
A

% Create the vector of loads
y = zeros(2*nn,1);
for k=1:nn
    % Enter the loads for node k into y
    y(2*k-1) = -loads(k,1);
    y(2*k)   = -loads(k,2);
end
y

forces = A\y