% problem 2.27
%
Aab = 1600 %mm
Aad = 3200
Ace = 3200
Abd = 2600
Aae = 3900
Abc = 1300
Agf = Aab
Adg = Aad
Acg = Aae
Aef = Abc
Adf = Abd
%
E = 200*10^8 %Pi
Lab = 4 %mM
Lad = 8
Lce = 8
Lbd = sqrt(4^2+8^2)
Lae = sqrt(12^2+6^2)
Lbc = sqrt(4^2+2^2)
Lgf = Lab
Ldg = Lad
Lcg = Lae
Lef = Lbc
Ldf = Lbd
%
ndof = 7*2
fdof = ndof-3
nelem = 11
Ksys = zeros(ndof)
%
% Set up vectors to assemble
Area = [Aab Agf Aad Adg Ace Abd Acg Aef Adf Aae Abc]./1000^2 %meters
Length = [Lab Lgf Lad Ldg Lce Lbd Lcg Lef Ldf Lae Lbc]
E = E*ones(1,length(Area))
phi = [pi/2 pi/2 0 0 0 pi-atan(4/8) pi-atan(6/12) pi-atan(2/4) atan(4/8) atan(6/12) atan(2/4)] 
Xinfo = [13 14 1 2;11 12 9 10;13 14 5 6;5 6 11 12;3 4 7 8;5 6 1 2;11 12 3 4;9 10 7 8;5 6 9 10;13 14 7 8; 1 2 3 4]
%%
Ksys = zeros(ndof)
for i=1:nelem
	k = swftrusstiff(Length(i),E(i),Area(i))
	T = swfgetT(phi(i),1)
	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 = zeros(fdof,1)
Pf(6) = -60*1000 %Newtons
%
% member forces
Pfm = zeros(fdof,1)
Psm = zeros(ndof-fdof,1)
%
% coordinate displacements
Df = inv(Kff)*(Pf - Kfs*Ds - Pfm)
%
% reaction forces
Ps = Ksf*Df + Kss*Ds + Psm
%
D = [Df' Ds']'
%D = [4.46 -.75 0 5.58 -.75 -13.74 -1.5 5.58 -5.96 -.75 -1.5 0 0 0]'./1000
% member forces
for i=1:nelem
	k = swftrusstiff(Length(i),E(i),Area(i))
	T = swfgetT(phi(i),1)
	Xmap = Xinfo(i,:)
	X = swfgetX(Xmap,ndof)
   X*D
   d = T*X*D
   s = k*d
   bars(i) = s(3)
	pause   
end
%
['The following forces are multiplied by P.']
['positive is tension']
['negative is compression']
bars ./1000