% Workshop Example 10  (WorkshopEx10.m)
%
% This script displays and analyzes a truss specified by a
% descriptor file containing the following information
% 
% 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.
%
% hooke - the vector of Hooke's (stiffness) constants for each 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.
%
% load - the matrix specifying the first set of loads on the
%      - nodes.  If there are n nodes, then this is an n-by-2
%      - matrix.  The entries in the first column are the 
%      - horizontal loads, and the entries in the second column
%      - are the vertical loads.
%
% 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 15, 2008

clear
clc

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]);
 
% Plot the bars
hold on
[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.25*xmargin, yn(k)+0.25*ymargin, num2str(k))
end
pause
 
% Set up the preliminary structural matrix
A0 = zeros(nb,2*nn);
for k=1:nb
    % Get the index for the first node (source node) of bar k
    
    % Get the index for the second node (target node) of bar k
    
    % Calculate the length of each bar
    
    % Calculate ck and sk
    
    
    % Enter the signed values of ck and sk into the appropriate
    % columns of row k.
   
    
    
    
end
A0
 
% Copy the columns that do not correspond to fixed nodes
index = 1;
for k=1:nn
    if (fixed(k,1) == 0)
        A(:,index) = A0(:,2*k-1);
        index = index+1;
    end
    if (fixed(k,2) == 0)
        A(:,index) = A0(:,2*k);
        index = index+1;
    end
end
% Copy the columns that do not correspond to fixed nodes
% and create the vector of external loads skipping over
% the loads that correspond to reactive loads
index = 1;
for k=1:nn
    if (fixed(k,1) == 0)
        A(:,index) = A0(:,2*k-1);
        f(index,1) = load(k,1);
        index = index+1;
    end
    if (fixed(k,2) == 0)
        A(:,index) = A0(:,2*k);
        f(index,1) = load(k,2);
        index = index+1;
    end
end
A
f

% Construct the stiffness matrix
C = diag(hooke);
K = A'*C*A
 
% Solve for the displacements that minimize the potential
% energy of the loaded truss using the LU factorization of
% the stiffness matrix
u = K \ f 
% Calculate the forces acting on the bars
fprintf('Force acting on each bar\n')
yb = C*A*u;

% Plot the bars color coded according to the forces
tload = sum(abs(load(:,1))) + sum(abs(load(:,2)))
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);
    if (abs(yb(k)) < 0.001*tload)
        plot([xn(sindx) xn(tindx)],[yn(sindx) yn(tindx)],'g')
    elseif (0 < yb(k))
        plot([xn(sindx) xn(tindx)],[yn(sindx) yn(tindx)],'b')
    else
        plot([xn(sindx) xn(tindx)],[yn(sindx) yn(tindx)],'r')
    end 
end
hold off

% Display the matrix of displacements for all the nodes
u0 = zeros(nn,2);
index = 1;
for k=1:nn
    if (fixed(k,1)==0)
        u0(k,1) = u(index);
        index = index+1;
    end
    if (fixed(k,2)==0)
        u0(k,2) = u(index);
        index = index+1;
    end
end
fprintf('Displacements for all nodes\n')
u0
 
% Calculate and display all the external forces including
% the reactive forces
f0 = A0'*yb;
ef = zeros(nn,2);
for k=1:nn
    ef(k,1) = f0(2*k-1);
    ef(k,2) = f0(2*k);
end
fprintf('External forces on all nodes\n')
ef