% problem 4.28
%
I = 260*10^-6 %m^4
% calculate twice -- once for each A below.
A = .03 %m^2
A = 100000000000000
E = 200*10^9 %Pa
L = [8 8 5 5 5 5]' %meters
type = 2 %frame
w = 18*1000 %N/m distributed load
%
ndof = 7*3
fdof = ndof-5
nelem = 6
%
% Set up vectors to assemble
phi = [0 0 -pi/2 -pi/2 -pi/2 -pi/2] 
Xinfo = [19 20 21 1 2 3;1 2 3 4 5 6;4 5 6 7 8 9;7 8 9 10 11 12;10 11 12 13 14 15;13 14 15 17 18 16]
%
Ksys = zeros(ndof)
for i=1:nelem
   k = swffrmstiff(L(i),E,I,A)
   pause
   T = swfgetT(phi(i),type)
   pause
	kg = T'*k*T
	Xmap = Xinfo(i,:)
	X = swfgetX(Xmap,ndof)
	Ksys = Ksys + X'*kg*X
	pause   
end
%%
%%
Kff = Ksys(1:fdof,1:fdof)
Kfs = Ksys(1:fdof,fdof+1:ndof)
Ksf = Ksys(fdof+1:ndof,1:fdof)
Kss = Ksys(fdof+1:ndof,fdof+1:ndof)
%
Ds = zeros(ndof-fdof,1)
%%%
%%%
% Assembly load vector
Pf = [0 0 0 0 0 -50 16 0 0 0 0 -30 0 0 0 0]'*1000
%
% member forces
% Assembly member load vector
Pmember = zeros(ndof,1)
Loads = [0 w*L(1)/2 w*L(1)^2/12 0 w*L(1)/2 -w*L(1)^2/12;
   0 0 0 0 0 0; 0 0 0 0 0 0;
   0 w*L(4)/2 w*L(4)^2/12 0 w*L(4)/2 -w*L(4)^2/12;
   0 0 0 0 0 0; 0 0 0 0 0 0] 
for i=1:nelem
	T = swfgetT(phi(i),type)
   Pm = T'*[Loads(i,:)]'
	Xmap = Xinfo(i,:)
	X = swfgetX(Xmap,ndof)
	Pmember = Pmember + X'*Pm
	pause   
end
%
%
Pfm = Pmember(1:fdof)
Psm = Pmember(fdof+1:ndof)
%
% coordinate displacements
Df = inv(Kff)*(Pf - Kfs*Ds - Pfm)
%
% reaction forces
Ps = Ksf*Df + Kss*Ds + Psm
%
D = [Df' Ds']'
% member forces
for i=1:nelem
	k = swffrmstiff(L(i),E,I,A)
	T = swfgetT(phi(i),type)
	Xmap = Xinfo(i,:)
	X = swfgetX(Xmap,ndof)
   X*D
   d = T*X*D
   s = k*d + Loads(i,:)'
   beams(i,:) = s'
	pause   
end
%
%

['The following displacements are in global coordinates.']
['Note that translations are close to zero.']
Df

['The following forces are in member coordinates. Column i = Member i.']
beams'./1000

['The reaction forces are in global coordinates']
Ps./1000