reads %truss2d.m data from an Excel file called 'truss2d_data.xls'
% Worksheets in the xls file have the same name same as the data read
clear all
% The output is written to the file output.txt
fid = fopen('output.m','w');
Nodes=xlsread('truss2d_data','Nodes');
Nelms=xlsread('truss2d_data','Nelms');
x=xlsread('truss2d_data','x');
NBC=xlsread('truss2d_data','NBC');
node=xlsread('truss2d_data','node');
props=xlsread('truss2d_data','props');
force=xlsread('truss2d_data','force');
% Print number of nodes and elements
fprintf(fid,'Number of nodes= %3i\n',Nodes);
fprintf(fid,'Number of elements= %3i\n',Nelms);
%Create array NN containg 1,2,3,..Nodes
NN=[1:1:Nodes];
dprint=[NN',x,NBC];
fprintf(fid,'Node numbers, their coordinates and BCs are:\n');
fprintf(fid,'%3i %12.8f %12.8f %3i %3i\n',dprint');
%Create array MM containg 1,2,3,..Nelms
MM=[1:1:Nelms];
%Augment [MM], [node] and [props] to form [pprint]
pprint=[MM',node,props];
fprintf(fid,'Element number, first and second nodes, A and E of elements are:\n');
fprintf(fid,'%3i %3i %3i %11.4e %11.4e\n',pprint');
lprint=[NN',force];
fprintf(fid,'Node number, x-force and y-force:\n');
fprintf(fid,'%3i %11.4e %11.4e\n',lprint');
% NEQ = Number of equations in [K] = No. of active DOFs
% First set NEQ=0
NEQ=0;
for i=1:Nodes
for j=1:2
if NBC(i,j)==1
NEQ=NEQ+1;
NBC(i,j)=NEQ;
end
end
end
fprintf(fid,'Number of equations= %3i\n',NEQ);
% Initialize [K] and {F}
K=zeros(NEQ);
F=zeros(NEQ,1);
% Form {F} from the matrix [force]
for i=1:Nodes
for j=1:2
ID=NBC(i,j);
if ID~=0
F(ID)=force(i,j);
end
end
end
% Compute [k] for each element and assemble [K]
for e=1:Nelms
for i=1:2
% Lv is the vector in the direction of the element
Lv(i)=x(node(e,2),i)-x(node(e,1),i);
end
% L(e) is the length of Element e, which is
%calculated as the norm of the vector {Lv}
L(e)=norm(Lv);
for i=1:2
% dc are the direction cosines
%dc(e,1) is the direction cosine l of Element e
%%dc(e,2) is the direction cosine m of Element e
dc(e,i)=Lv(i)/L(e);
%call end
function stiff1.m to compute [k] for Element e
[k]=stiff1(L(e),dc(e,:),props(e,:));
for m=1:2
for n=1:2
%scan each row of [k]
lrow=(m-1)*2+n;
% lrow is the local row number in [k]
grow=NBC(node(e,m),n);
% grow is the corresponding global row
if grow~=0
% if grow=0 strike the row, do not assemble
for i=1:2
for j=1:2
lcol=(i-1)*2+j;
gcol=NBC(node(e,i),j);
if gcol~=0
% if gcol=0 strike the column
K(grow,gcol)=K(grow,gcol)+k(lrow,lcol);
end
end
end
end
end
end
end
% solve for{Q}
Q=inv(K)*F;
% Initialize the array u(i,j) of size Nodes x 2
u=zeros(Nodes,2);
for i=1:Nodes
for j=1:2
if NBC(i,j)~=0
u(i,j)=Q(NBC(i,j));
end
end
end
for e=1:Nelms
for i=1:2
del(i)=u(node(e,2),i)-u(node(e,1),i);
end
%del is a vector of relative displacements
DeltaL=sum(dc(e,:).*del);
% DeltaL is computed as the dot product of
% del and direction cosines vector
P(e,1)=props(e,1)*props(e,2)*DeltaL/L(e);
end
%Create array NN containg 1,2,3,..Nodes
NN=[1:1:Nodes];
%Create array MM containing 1,2,3, ...Nelms
MM=[1:1:Nelms];
%Augment [NN] and [u] to form [uprint]
uprint=[NN',u];
%Augment [MM] and [P] to form [Pprint]
Pprint=[MM',P];
% print displacements[uprint]
fprintf(fid,'Nodal displacements:\n');
fprintf(fid,'%3i %11.4e %11.4e\n',uprint');
% Print element forces [Pprint]
fprintf(fid,'Element forces:\n');
fprintf(fid,'%3i %11.4e\n',Pprint');
fclose(fid);
mf=input('Enter a Magnification factor, usually close 10')
zfactor=mf*min(L)/100/(norm(Q)/length(Q));
for e=1:Nelms
for i=1:2
xp(i)=x(node(e,i),1);
yp(i)=x(node(e,i),2);
end
plot(xp,yp)
hold on
end
for e=1:Nelms
for i=1:2
xp(i)=x(node(e,i),1)+u(node(e,i),1)*zfactor;
yp(i)=x(node(e,i),2)+u(node(e,i),2)*zfactor;
end
plot(xp,yp,'r')
hold on
end