%=====CME434 Fall 06: Solution to HW4 by TSD 10/13/06
%Begin with your script from HW3.
%10 (Deferred from HW3) Write a script to open the input data file hw4data.txt (the same, for now, as hw3data.txt) and create an initial workspace like the one opened at the start of HW3.  Put this script at the beginning of your HW3 script, instead of load(). Now you will start each problerm with a cleared workspace and create the initial workspace from a problem data file. Test your script to see that the initial workspace is correct.
status = fclose('all')		%Make sure all files are closed
%File names are hw4data.txt hw41data.txt & hw42data.txt
FIDI_NAME = 'hw4data.txt'	%Name of FIDI

%----- Read Input Data File into Row Vectors A & B -----
FIDI = fopen(FIDI_NAME, 'r');	%Open input file
disp(blanks(1));		%Write blanks to screen
TEXT = fgets(FIDI)		%Read captions, 4 lines
TEXT = fgets(FIDI)
TEXT = fgets(FIDI)
TEXT = fgets(FIDI)
DATA = str2num(fgets(FIDI));
[nn ne nm ndim nep nqn] = deal(DATA(1), DATA(2), DATA(3), DATA(4), DATA(5), DATA(6))
TEXT = fgets(FIDI)
DATA = str2num(fgets(FIDI));
[nq nf nmpc] = deal(DATA(1), DATA(2), DATA(3));

TEXT = fgets(FIDI)
XYZ = zeros(nn,3);		%make node coordinates table
for i=1:nn
 DATA = str2num(fgets(FIDI));
 for j=1:ndim
  XYZ(DATA(1),j) = DATA(1+j);
 end
end

TEXT = fgets(FIDI)
CTAB = zeros(ne,2);		%make CTAB: connectivity table
AE = zeros(ne,1); 		%make AE: members property table
for i=1:ne
 DATA = str2num(fgets(FIDI));
 [CTAB(DATA(1),1) CTAB(DATA(1),2) AE(DATA(1))] = deal(DATA(2), DATA(3), DATA(4));
end

TEXT = fgets(FIDI)
QTAB = zeros(nq,2);		%make QTAB: essential displacements table
for i=1:nq
 DATA = str2num(fgets(FIDI));
 [QTAB(i,1) QTAB(i,2)] = deal(DATA(1), DATA(2));
end

TEXT = fgets(FIDI)
FTAB = zeros(nf,2);		%make FTAB: specified forces table
for i=1:nf
 DATA = str2num(fgets(FIDI));
 [FTAB(i,1) FTAB(i,2)] = deal(DATA(1), DATA(2));
end

TEXT = fgets(FIDI)
MATERIAL = zeros(nm,2);		%make FTAB: specified forces table
for i=1:nm
 DATA = str2num(fgets(FIDI));
 nmp=1;				%number of material properties
 for j=1:nmp
  MATERIAL(DATA(1),j) = DATA(1+j);
 end
end

TEXT = fgets(FIDI)
status = fclose('all');

%1
ns = nn*nqn

%2  Write a script to calculate bar lengths Lx, Ly, L, l=cos(theta) and m=sin(theta) for each of the ne elements.
%Store the results in (1 x ne) row vectors [L.e], [l.e] and [m.e].
for e=1:ne
    Lx=XYZ(CTAB(e,2),1)-XYZ(CTAB(e,1),1)
    Ly=XYZ(CTAB(e,2),2)-XYZ(CTAB(e,1),2)
    L(e)=sqrt(Lx*Lx+Ly*Ly)
    l(e)=Lx/L(e)
    m(e)=Ly/L(e)
end

%3  Write a script to draw the truss by plotting each element (graph the endpoints) in the same figure.
figure(1);
axis equal;
 Xmin=XYZ(CTAB(1,1),1);
 Ymin=XYZ(CTAB(1,1),2);
 Xmax=Xmin;
 Ymax=Ymin;
hold on
for e=1:ne
    Xe2=XYZ(CTAB(e,2),1);
    Xe1=XYZ(CTAB(e,1),1);
    Ye2=XYZ(CTAB(e,2),2);
    Ye1=XYZ(CTAB(e,1),2);
    plot([Xe1 Xe2], [Ye1 Ye2]);
    Xmin=min([Xmin Xe1 Xe2]);
    Xmax=max([Xmax Xe1 Xe2]);
    Ymin=min([Ymin Ye1 Ye2]);
    Ymax=max([Ymax Ye1 Ye2]); 
end
e=max(Xmax-Xmin,Ymax-Ymin)/32;	%add border space
plot(Xmin-e,Ymin-e);
plot(Xmax+e,Ymax+e);

%Label nodes
for i=1:nn
    text(XYZ(i,1), XYZ(i,2), num2str(i));
end
hold off

%4  Write a script to calculate k=AE/L for each element.
%Store the results in (l x ne) row vector [k.e]
for e=1:ne
    k(e)=AE(e)/L(e);
end
k

%5  Write a script, given nf and FTAB, to initialize the column vector F with the specified applied forces.
%Note: an applied force is zero unless some other value is specified.
F=zeros(ns,1);
for i=1:nf
    F(FTAB(i,1))=FTAB(i,2);
end
F

%6  Write a Script, given nq and QTAB, to initialize the column vector Q with the specified displacements.
%Note: To create Q a value must be given. For now, let the Q(i) be zero except where another value is given.
Q=zeros(ns,1);
for i=1:nq
    Q(QTAB(i,1))=QTAB(i,2);
end
Q

%7  Write a script to calculate Ke and Ge for each element.
%Assemble the results in K using K=Sum(Ge'*Ke*Ge)
K=zeros(ns,ns);
Ge=zeros(2*nqn,ns);
for e=1:ne
    M=[l(e)*l(e) l(e)*m(e);l(e)*m(e) m(e)*m(e)]
    Ke=k(e)*[M -M;-M M]
    Ge=0*Ge;
    Ge(1,2*CTAB(e,1)-1)=1;
    Ge(2,2*CTAB(e,1))=1;
    Ge(3,2*CTAB(e,2)-1)=1;
    Ge(4,2*CTAB(e,2))=1
    K=K+Ge'*Ke*Ge
end

%8  Write a script to apply the boundary conditions and solves for displacements.
%Assume that Q(1) ... Q(nq) are the known displacements. Some may be non-zero.
%Partition Q, K, and F, solve for the free displacements QF=inv(KF)*(FF-KEF'*QE.
QE=Q(1:nq)
%QF=Q(nq+1:ns)
KE=K(1:nq,1:nq)
KF=K(nq+1:ns,nq+1:ns)
KEF=K(1:nq,nq+1:ns)
%FE=F(1:nq)
FF=F(nq+1:ns)
QF=inv(KF)*(FF-KEF'*QE)

%9  Write a Script to calculate the Reactions FE=KE*QE+KEF*QF.
%Complete Q and F by including the calculated values of QF and FE.
FE=KE*QE+KEF*QF
Q(nq+1:ns)=QF
F(1:nq)=FE

%11  Write a script to calculate the nodal forces for each element Fe=Ke*Qe 
%Save the bar force for each element in col. vector P. P(e)=sqrt(Fe2x^2+Fe2y^2)
P=zeros(ne,1)
for e=1:ne
    M=[l(e)*l(e) l(e)*m(e);l(e)*m(e) m(e)*m(e)];
    Ke=k(e)*[M -M;-M M]
    Ge=0*Ge;
    Ge(1,2*CTAB(e,1)-1)=1;
    Ge(2,2*CTAB(e,1))=1;
    Ge(3,2*CTAB(e,2)-1)=1;
    Ge(4,2*CTAB(e,2))=1
    Fe=Ke*Ge*Q
    P(e)=[l(e) m(e)]*Fe(3:4)
end

%11  Write a script to calculate the nodal force for each element Fe=Ke*Qe 
%Save the bar force for each element in col. vector P. P(e)=sqrt(Fe2x^2+Fe2y^2)
P=zeros(ne,1)
for e=1:ne
    M=[l(e)*l(e) l(e)*m(e);l(e)*m(e) m(e)*m(e)];
    Ke=k(e)*[M -M;-M M]
    Ge=0*Ge;
    Ge(1,2*CTAB(e,1)-1)=1;
    Ge(2,2*CTAB(e,1))=1;
    Ge(3,2*CTAB(e,2)-1)=1;
    Ge(4,2*CTAB(e,2))=1
    Fe=Ke*Ge*Q
    P(e)=[l(e) m(e)]*Fe(3:4)
end


%12 Rework your part %3 graphing script to also display the node numbers and the deflected shape of the truss.  Draw the truss as before. Use text() to show the node number on the graph at the coordinates of each node. Then superimpose the deformed truss shown in dash lines. For this, simply plot again, with displacements added to the original coordinates to locate the displaced coordinates.  To be distinct, the displacements will need to be maginfied by a factor of about 32 (experiment a little to see what you like) before adding them to the coordinates.  Graduate students will also draw the deformed truss with tension members in green, compression members in red.
figure(2);
axis equal;
%draw original truss
 Xmin=XYZ(CTAB(1,1),1);
 Ymin=XYZ(CTAB(1,1),2);
 Xmax=Xmin;
 Ymax=Ymin;
hold on
for e=1:ne
    Xe2=XYZ(CTAB(e,2),1);
    Xe1=XYZ(CTAB(e,1),1);
    Ye2=XYZ(CTAB(e,2),2);
    Ye1=XYZ(CTAB(e,1),2);
    plot([Xe1 Xe2], [Ye1 Ye2]);
    Xmin=min([Xmin Xe1 Xe2]);
    Xmax=max([Xmax Xe1 Xe2]);
    Ymin=min([Ymin Ye1 Ye2]);
    Ymax=max([Ymax Ye1 Ye2]); 
end
e=max(Xmax-Xmin,Ymax-Ymin)/32;	%add border space
plot(Xmin-e,Ymin-e);
plot(Xmax+e,Ymax+e);

%Label nodes
for i=1:nn
    text(XYZ(i,1), XYZ(i,2), num2str(i));
end

%draw deformed truss
mag=32
for e=1:ne
    n1=CTAB(e,1);
    n2=CTAB(e,2);
    Xe2=XYZ(n2,1)+mag*Q(2*n2-1);
    Xe1=XYZ(n1,1)+mag*Q(2*n1-1);
    Ye2=XYZ(n2,2)+mag*Q(2*n2);
    Ye1=XYZ(n1,2)+mag*Q(2*n1);
    if P(e)<0
        plot([Xe1 Xe2], [Ye1 Ye2], 'r--');  %red compression bars
    else
        plot([Xe1 Xe2], [Ye1 Ye2], 'g--');  %green tension bars
    end
end
hold off

%13 In general, hw4data.txt is a problem definition file for 2D trusses used to present data to your script. The present version of hw4data.txt was copied from hw3data.txt.  Delete the current hw4data.txt, copy hw41data.txt and rename the copy hw4data.txt.  Run your script on this new problem (chevron bracing in a building core to resist wind loads).  Do the same with hw42data.txt (a bridge with a truck at midspan, landscape mode where y is up).  Your script should solve all three problems without modification.  If it does not, modify your script (most likely by substituting variables for literals) until the same script works on all three problems.