clear; close all
%
%  Problem 3.24 thru 3.26
%
%  Let A = 1 in^2
%  E = 10 million psi
%
%
%  Step 1: define nodes & elems
%
nodes = 

nnodes =
ndof = 2*nnodes;

elems =

nelems =

%
%  Step 2: element properties
%
A = 1;
E = 10.e6;
%
%  Step 3: loads
%
loads = [2  2  15000     %% node 2 y dir
         3  2  15000];   %% node 3 y dir
nloads = size(loads,1);
%
%  Step 4: Assembly
%
bigk = zeros(ndof,ndof);
for i=1:nelems
    n1 =
    n2 =
    k = truss2d(n1,n2,nodes,A,E);
    dof = [2*n1-1 2*n1 2*n2-1 2*n2];
  bigk(dof,dof) = bigk(dof,dof) + k;
end;
%
%  Step 6: global load vector
%
bigp = zeros(ndof,1);
dof = 2*loads(:,1) -2 + loads(:,2);
bigp(dof) = loads(:,3);
%
%  Step 7: specify supports and solve
%
%
bc = [1  1   %% node 1, dx
      1  2   %% node 1, dy
      4  1   %% node 4, dx
      4  2]; %% node 4, dy
bcdof = 2*bc(:,1) -2 + bcdof(:,2);
freedof = 1:ndof;
freedof(bcdof) = [];
bigd = zeros(ndof,1);
bigd(freedof) = bigk(freedof,freedof)\bigp(freedof);
%
%  Step 8: tabulate global displ
%
fprintf('Displacements (mm):\n');
for i=1:nnodes
    dof = [2*i-1 2*i];
    fprintf('%4d %8.1f  %8.1f\n',i,bigd(dof)*1000);
end
%
%  Step 9: reactions and loads
%
bigf = bigk*bigd;
fprintf('Reactions & loads (kN):\n');
for i=1:nnodes
    dof = [2*i-1 2*i];
    fprintf('%4d %8.3f  %8.3f\n',i,bigf(dof)/1000);
end
%
%  Step 10: element results (later)
%
fprintf('Element stresses: (MPa)\n');
for i=1:nelems
    sig = truss2s(nodes,n1,n2,A,E,bigd);
    fprintf('%3d %8.3\n',sig/1.e6);
end