Truss Program Listing:

ENGN1300: Mechanics of Structures
FEM Truss Analysis Program Listing.

Graphics Version

> restart:

# Finte Element Truss Analysis Program

#********** GIVE THE NAME OF THE INPUTFILE IN THE LINE BELOW ***********

# Use frontslash/ in place of backslash\

# Enclose the name in single backquotes`

# End the line with a colon:

infilename:=

`C:/Documents and Settings/laptop306/Desktop/tutorial.inp.txt`:

#*************************** THEN HIT RETURN ***************************

#***********************************************************************

# load plotting and linear algebra maple packages

with(plots):with(linalg):with(plottools):

# define an output file name

lifn:=length(infilename):

outfilename:=cat(substring(infilename,1..lifn-8),`.out.txt`):

# ================= Read data from the input file ==================

infile:=fopen(infilename,READ):

# Title:

runtitle := readline(infile):

# No. nodes and nodal coordinates

nnode := fscanf(infile,`%d`)[1]:

coord := array(1..nnode,1..2):

for i from 1 to nnode do

coord[i,1] := fscanf(infile,`%f`)[1]:

coord[i,2] := fscanf(infile,`%f`)[1]:

od:

# No. elements connectivity, and material data.

# connect[i,1] and connect[i,2] hold the two joint numbers for the

# ith member. connect[i,3] holds EA for that member

nelem := fscanf(infile,`%d`)[1]:

connect := array(1..nelem,1..3):

for i from 1 to nelem do

for j from 1 to 2 do

connect[i,j] := fscanf(infile,`%d`)[1]:

od:

connect[i,3] := fscanf(infile,`%e`)[1]:

od:

# No. degrees of freedom with prescribed displacements; prescribed

# displacements

nfix := fscanf(infile,`%d`)[1]:

fixnodes := array(1..nfix,1..3):

for i from 1 to nfix do

# fixnodes[i,1] holds the number of the ith fixed node

# fixnodes[i,2] holds the direction (1 or 2)

# fixnodes[i,3] holds the valus of the prescribed displacement

fixnodes[i,1] := fscanf(infile,`%d`)[1]:

fixnodes[i,2] := fscanf(infile,`%d`)[1]:

fixnodes[i,3] := fscanf(infile,`%e`)[1]:

od:

# No. loaded nodes, with the loads

nload := fscanf(infile,`%d`)[1]:

loads := array(1..nload,1..3):

# Node number of ith loaded node is in loads[i,1]

# F1 and F2 for the ith loaded node are in loads[i,2] and loads[i,3]

for i from 1 to nload do

loads[i,1] := fscanf(infile,`%d`)[1]:

loads[i,2] := fscanf(infile,`%e`)[1]:

loads[i,3] := fscanf(infile,`%e`)[1]:

od:

fclose(infile):

#********************************GRAPHICS************************************************

# Plot the undeformed mesh

elplot:=[seq([[coord[connect[i,1],1],coord[connect[i,1],2]],

[coord[connect[i,2],1],coord[connect[i,2],2]]], i=1..nelem)]:

pnodes:=textplot({seq([coord[nn,1],coord[nn,2],convert(nn,string)],nn=1..nnode)},

font=[TIMES,BOLD,16]):

pels:=textplot({seq([(coord[connect[ne,1],1]+coord[connect[ne,2],1])/2,

(coord[connect[ne,1],2]+coord[connect[ne,2],2])/2,convert(ne,string)],

ne=1..nelem)},font=[HELVETICA,12]):

pm:=plot(elplot,axesfont=[TIMES,ROMAN,6],

color=magenta,labelfont=[TIMES,ROMAN,12], titlefont=[TIMES,ROMAN,14],title=cat(runtitle,`: Mesh`),

thickness=3,labels=[`x1`,`x2`],scaling=constrained):

display(pm,pnodes,pels);

# Plot loads and BC

# Scale the load vector by max element length*max load magnitude*.2

maxlength:= max(seq(sqrt((coord[connect[ie,1],1]-coord[connect[ie,2],1])^2+

(coord[connect[ie,1],2]-coord[connect[ie,2],2])^2),ie=1..nelem)):

maxf:=max(seq(sqrt(loads[ii,2]^2+loads[ii,3]^2),ii=1..nload)):

loadscale:=0.2*maxlength/(maxf+0.00001*maxlength):

# graphics of laod vectors

loadarrows:=seq(arrow([coord[loads[il,1],1],coord[loads[il,1],2]], [coord[loads[il,1],1]+loadscale*loads[il,2],coord[loads[il,1],2]+loadscale*loads[il,3]],

.02, .1, .3, color=green),il=1..nload):

# Scale the displacement arrows by max element length*max presc. displacement*.2

maxd:=max(seq(abs(fixnodes[ii,3]),ii=1..nfix)):

bcscale:=0.2*maxlength/(maxd+0.00001*maxlength):

# graphics of BC vectors

bcarrows:=seq(arrow(

[coord[fixnodes[il,1],1]-(2-fixnodes[il,2])*(fixnodes[il,3]+.05*maxlength),

coord[fixnodes[il,1],2]-(fixnodes[il,2]-1)*(fixnodes[il,3]+.05*maxlength)],

[coord[fixnodes[il,1],1],coord[fixnodes[il,1],2]],.01,.03,1-bcscale*fixnodes[il,3], color=blue),il=1..nfix):

# graphics of load labels

plls:=textplot({seq([(coord[loads[il,1],1]+loadscale*loads[il,2]/2,coord[loads[il,1],2]+loadscale*loads[il,3]/2,

substring(convert(sqrt(loads[il,2]^2+loads[il,3]^2),string),1..4))],

il=1..nload)},font=[HELVETICA,12]):

display(plls,pm,bcarrows,loadarrows,title=cat(runtitle,`: Load and BC`));

#********************************END OF GRAPHICS***********************************

Warning, the name changecoords has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name arrow has been redefined

> #

#================ ELEMENT STIFFNESS MATRIX ===================

elstif := proc (xa,ya,xb,yb,EA)

local elstif, nx, ny, length;

# Define the element length and direction vector

length := sqrt((xa-xb)*(xa-xb)+(ya-yb)*(ya-yb));

nx := (xb-xa)/length;

ny := (yb-ya)/length;

# Define the element stiffness

elstif := evalm((EA/length)*array(1..4,1..4,

[[nx*nx,nx*ny,-nx*nx,-nx*ny],

[ny*nx,ny*ny,-ny*nx,-ny*ny],

[-nx*nx,-nx*ny,nx*nx,nx*ny],

[-ny*nx,-ny*ny,ny*nx,ny*ny]]));

end:

#===========Assemble the global stiffness matrix==============

Stif := array(1..2*nnode,1..2*nnode):

for i from 1 to 2*nnode do

for j from 1 to 2*nnode do

Stif[i,j] := 0.0:

od:

# Loop over all the elements

for lmn from 1 to nelem do

# Set up the stiffness for the current element

a := connect[lmn,1]:

b := connect[lmn,2]:

EA := connect[lmn,3]:

k := array(1..4,1..4):

k:=elstif(coord[a,1],coord[a,2],coord[b,1],coord[b,2],EA):

# Add the current element stiffness to the global stiffness

for i from 1 to 2 do

for ii from 1 to 2 do

for j from 1 to 2 do

for jj from 1 to 2 do

rw := 2*(connect[lmn,i]-1)+ii:

cl := 2*(connect[lmn,j]-1)+jj:

Stif[rw,cl] := Stif[rw,cl] + k[2*(i-1)+ii,2*(j-1)+jj]:

od:

# ==================== Assemble global residual vector ============

# Define the residual

resid := array(1..2*nnode):

for i from 1 to 2*nnode do

resid[i] := 0.:

od:

# Loop loaded nodes

for i from 1 to nload do

a := loads[i,1]:

f1 := loads[i,2]:

f2 := loads[i,3]:

resid[2*a-1] := resid[2*a-1] + f1:

resid[2*a] := resid[2*a] + f2:

od:

# Modify the global stiffness and residual to include constraints

for i from 1 to nfix do

rw := 2*fixnodes[i,1]+fixnodes[i,2]-2:

for j from 1 to 2*nnode do

resid[j] := resid[j] - Stif[j,rw]*fixnodes[i,3]:

Stif[rw,j] := 0:

Stif[j,rw] := 0:

od;

Stif[rw,rw] := 1.0:

resid[rw] := fixnodes[i,3]:

od:

# ================== Solve the FEM equations ===================

u := linsolve(Stif,resid):

# Define a procedure to calculate the element strains

elstrain := proc (xa,ya,xb,yb,uax,uay,ubx,uby)

local elstrain, nx, ny, length, uel;

# length and direction vector for the element

length := sqrt((xa-xb)*(xa-xb)+(ya-yb)*(ya-yb));

nx := (xb-xa)/length;

ny := (yb-ya)/length;

# Element displacement vector

uel := array(1..4,[uax,uay,ubx,uby]);

# Element strains

elstrain := ((ubx-uax)*nx+(uby-uay)*ny)/length;

end:

# compute element forces and strains:

for lmn from 1 to nelem do

a := connect[lmn,1]:

b := connect[lmn,2]:

EA := connect[lmn,3]:

xa := coord[a,1]:

ya := coord[a,2]:

xb := coord[b,1]:

yb := coord[b,2]:

uxa := u[2*a-1]:

uya := u[2*a]:

uxb := u[2*b-1]:

uyb := u[2*b]:

Strain[lmn] := elstrain(xa,ya,xb,yb,uxa,uya,uxb,uyb):

Force[lmn]:=EA*Strain[lmn];

od:

#********************************GRAPHICS*****************************************

#=====Plot the deformed mesh===

#=====Displacements are scaled by 0.1*maximum displacement*maximum element length

maxu:=max(seq(abs(u[ii]),ii=1..2*nnode)):

plotscale:=0.1*maxlength/maxu:

print(`Displacements scaled by:`,plotscale);

dcoord:=array(1..nnode,1..2):

for a from 1 to nnode do

xa := coord[a,1]:

ya := coord[a,2]:

uxa := u[2*a-1]:

uya := u[2*a]:

dcoord[a,1]:=xa+plotscale*uxa:

dcoord[a,2]:=ya+plotscale*uya:

od:

defplot:=[seq([[dcoord[connect[i,1],1],dcoord[connect[i,1],2]],

[dcoord[connect[i,2],1],dcoord[connect[i,2],2]]], i=1..nelem)]:

pd:=plot(defplot,axesfont=[TIMES,ROMAN,10],

color=green,labelfont=[TIMES,ROMAN,14], titlefont=[TIMES,BOLD,16],title=`Truss Mesh`,

thickness=3,labels=[`x1`,`x2`],scaling=constrained,title=cat(runtitle, `: Deformation`)):

dloadarrows:=seq(arrow([dcoord[loads[il,1],1],dcoord[loads[il,1],2]],

[dcoord[loads[il,1],1]+loadscale*loads[il,2],dcoord[loads[il,1],2]+loadscale*loads[il,3]],

.02, .1, .3, color=green),il=1..nload):

# Scale the bc vector by max element length*max presc. displacement*.2

maxd:=max(seq(abs(fixnodes[ii,3]),ii=1..nfix)):

bcscale:=0.2*maxlength/(maxd+0.00001*maxlength):

# graphics of bc

dbcarrows:=seq(arrow(

[dcoord[fixnodes[il,1],1]-(2-fixnodes[il,2])*(fixnodes[il,3]+.05*maxlength),

dcoord[fixnodes[il,1],2]-(fixnodes[il,2]-1)*(fixnodes[il,3]+.05*maxlength)],

[dcoord[fixnodes[il,1],1],dcoord[fixnodes[il,1],2]],.01,.03,1-bcscale*fixnodes[il,3], color=blue),il=1..nfix):

display(pd,pm,dbcarrows,dloadarrows);

#====Display element forces graphically

pforce:=textplot({seq([(coord[connect[ne,1],1]+coord[connect[ne,2],1])/2,(coord[connect[ne,1],2]+coord[connect[ne,2],2])/2, convert(evalf(Force[ne],3),string)],ne=1..nelem)},font=[HELVETICA,14]):

display({pm,pnodes,pforce},titlefont=[TIMES,BOLD,16],title=cat(runtitle,`: Forces`),

thickness=3,scaling=constrained,axes=none);

#******************************** END OF GRAPHICS ***********************************

# Print nodal displacements, element strains and forces

outfile:=fopen(outfilename,WRITE):

fprintf(outfile, `%s\n`,runtitle):

fprintf(outfile,`%s %3d %3d\n` ,`Number of Nodes and Elements:`,nnode, nelem):

fprintf(outfile,`%s\n`,`Node# x1 x2`):

for i from 1 to nnode do

fprintf(outfile,`%5d %9.4f %9.4f\n`,

i,coord[i,1],coord[i,2]):

od:

fprintf(outfile,`%s\n`,`El# Node 1 Node 2 EA`):

for i from 1 to nelem do

fprintf(outfile,`%4d %10d %10d %9.4f\n`,

i,connect[i,1],connect[i,2],connect[i,3]):

od:

fprintf(outfile,`%s\n`,`Fixed Node# DOF Value`):

for i from 1 to nfix do

fprintf(outfile,`%9d %12d %10.4f\n`,

fixnodes[i,1],fixnodes[i,2],fixnodes[i,3]):

od:

fprintf(outfile,`%s\n`,`Loaded Node# P1 P2`):

for i from 1 to nload do

fprintf(outfile,`%9d %13.4f %9.4f\n`,

loads[i,1],loads[i,2],loads[i,3]):

od:

fprintf(outfile,`%s\n`,` `):

fprintf(outfile,`%s\n`,`*********************************`):

fprintf(outfile,`%s\n`,` `):

fprintf(outfile,`%s\n`,`Nodal Displacements:`):

fprintf(outfile,`%s\n`,` Node u1 u2`):

for i from 1 to nnode do

fprintf(outfile,`%3d %8.4f %8.4f\n`,i,u[2*i-1],u[2*i]):

od:

fprintf(outfile,`%s\n`,`Strains and Forces:`):

fprintf(outfile,`%s\n`,`Element Strain Force`):

for i from 1 to nelem do

fprintf(outfile,`%3d %9.4f %9.4f\n`,i,Strain[i],Force[i]):

od:

fclose(outfile):

print(`End`);