TRUSS ANALYSIS v4.1 was first created in 14 June 2001 using QBasic and updated
in 09 December 2023 using Phyton by Ernel Solis Filipinas
(ernel_filipinas@yahoo.com).
https://colab.research.google.com/drive/1xNp5Hd96YDruSFZOiwcobsjLhk-DI_Aa?usp=sharing
(For educational
reference only)
-----------------------------------------------------copy
below-----------------------------------------------------
import math
import numpy as np
## INPUT DATA
#From Cañete (2021)
Example 3-3 of pagse 98-104.
# Define joint
coordinate [x,y] in mm
jdat1 = np.array([
[0,0],
[3900,3500],
[4800,0],
[9000,0]])
# Define node/joint
type [] where 1-Hinge 2-Roller x 3-Roller y 4-Free Joint
jdat3 =
np.array([[1],[4],[1],[1]])
# Define horizontal
load in kN
jdat4 =
np.array([[0],[42],[0],[0]])
# Define vertical
load in kN
jdat5 =
np.array([[0],[-56],[0],[0]])
# Define
element/member incidences
mdat1 =
np.array([[1,2],[2,3],[4,2]])
# Define Area in
mm2
mdat3 =
np.array([[2500],[1000],[4500]])
# Define Modulus of
Elasticity (GPa)
mdat4 =
np.array([[200],[200],[200]])
#----------------------------------------------------------------------------------
#----------------------------------------------------------------------------------
## CALCULATION
FUNCTION
# Step
0---------------------------------------------------------------------------
def
title_disclaimer():
print("TRUSS ANALYSIS v4.1 first created in 14 June 2001 using
QBasic")
print("Updated in 09 December 2023 using Phyton")
print("by Ernel Solis Filipinas. \ne:
ernel_filipinas@yahoo.com\n\n")
print("Disclaimer:\n")
print("No Liability is accepted by its software authors for any direct,
")
print("indirect, consequential or incidental loss or damage arising out of
")
print("the software use or any mistakes and negligence in developing this
software.")
print("The organisation or person using the software bears all risks and
")
print("responsibility for the quality and performance of the software.
\n")
print("Reference:\n")
print("Espejo, Isagani S. (2000). The Matrix Method - Computer-Aided
Structural Analysis")
print("Cañete, Alberto (2021). Principles of Matrix Structural
Analysis\n")
#Step
1---------------------------------------------------------------------------
# Define functions
checking number of joint and number of member
def
check(jdat1,jdat3,jdat4,jdat5,mdat1, mdat3, mdat4):
print("Step 1: Checking number of joint and number of
member-----------------------\n")
noj =
jdat1.shape[0]
nojdat3 = jdat3.shape[0]
nojdat4 = jdat4.shape[0]
nojdat5 = jdat5.shape[0]
print
("Number of joints =" ,noj)
if
noj == nojdat3:
if noj == nojdat4:
if noj == nojdat5:
print("The number joints and
joint datum is OK")
else:
print("Number of joints not
equal to number of vertical loads")
else:
print("Number of joints not equal to number of
horizontal loads ")
else:
print("Number of joints not equal to number of joint
type")
nom =
mdat1.shape[0]
nonmdat3 = mdat3.shape[0]
nonmdat4 = mdat4.shape[0]
print
("Number of members " ,nom)
if
nom == nonmdat3:
if nom == nonmdat4:
print("The number members and member datum is
OK\n")
else:
print("Number of members not equal to number
of member areas\n")
else:
print("Number of members not equal to number of member
modulus\n")
#Step
2---------------------------------------------------------------------------
# Define functions
for NDF, TDOF
def
DOF(jdat1,jdat3):
print("Step 2: Compute for NDF & TDOF-----------------------\n")
TDOF
= 0
noj =
jdat1.shape[0]
NDF =
np.zeros(shape=(noj, 2),dtype=int)
for i
in range (0,noj):
print("Joint: ", i+1)
if jdat3[i] == 1:
NDF[i,0] = 0
NDF[i,1] = 0
elif jdat3[i] == 2:
NDF[i,0] = TDOF + 1
NDF[i,1] = 0
TDOF+=1
elif jdat3[i] == 3:
NDF[i,0] = 0
NDF[i,1] = TDOF + 1
TDOF+=1
elif jdat3[i] == 4:
NDF[i,0] = TDOF + 1
NDF[i,1] = TDOF + 2
TDOF+=2
else:
print("Please check the joint data")
print("NDF Array x = ", NDF[i,0],"NDF Array y =
", NDF[i,1])
print
("TDOF = ",TDOF,"\n")
return NDF, TDOF
#Step
3---------------------------------------------------------------------------
# Define functions
for Qsys
def
Q(jdat1,jdat4,jdat5, NDF,TDOF):
print("Step 3: Compute for Qsys-----------------------\n")
noj =
jdat1.shape[0]
Qsys
= np.zeros(shape=(TDOF),dtype = int)
for i
in range (0,noj):
for j in range (0,2):
if NDF[i,j] > 0:
if j==0:
Qsys[NDF[i,0]-1] =
jdat4[i]
print("Qsys",NDF[i,0],"=", Qsys[NDF[i,0]-1],"kN")
if j==1:
Qsys[NDF[i,1]-1] =
jdat5[i]
print("Qsys",NDF[i,1],"=", Qsys[NDF[i,1]-1],"kN")
print("")
return Qsys
#Step
4---------------------------------------------------------------------------
# Define functions
for Mcode
def MCODE(mdat1):
print("Step 4: Compute for Mcode-----------------------\n")
nom =
mdat1.shape[0]
mcode
= np.zeros(shape=(nom,4),dtype = int)
MDIFF
= 0
for i
in range(0,nom):
print("Member:",i+1)
print("Member Incidence",mdat1[i,0],mdat1[i,1])
mcode[i,0] = NDF[mdat1[i,0]-1,0]
mcode[i,1] = NDF[mdat1[i,0]-1,1]
mcode[i,2] = NDF[mdat1[i,1]-1,0]
mcode[i,3] = NDF[mdat1[i,1]-1,1]
print("mcode",mcode[i,0],mcode[i,1],mcode[i,2],mcode[i,3])
MAX = 0
MIN = TDOF
for j in range (0,4):
if mcode[i,j] == 0:
continue
elif mcode[i,j] > MAX:
MAX = mcode[i,j]
if mcode[i,j] < MIN:
MIN = mcode[i,j]
print("MAX",MAX)
print("MIN",MIN)
if (MAX-MIN) >= MDIFF:
MDIFF = MAX - MIN
print("MDIFF",MDIFF,"\n")
HBW =
MDIFF+1
print("HBW",HBW,"\n")
return mcode, HBW
#Step
5---------------------------------------------------------------------------
# Define functions
for Ksys
def
create_stiffness_matrix(jdat1, mdat1, mdat3, mdat4,TDOF,mcode):
print("Step 5: Compute for Ksys-----------------------\n")
#
Initialize matrix
nom =
mdat1.shape[0]
K =
np.zeros((4,4))
Ksys
= np.zeros((TDOF,TDOF))
mprop
= np.zeros((nom,4))
#
Extract node coordinates=
#x=jdat1[:,0] # sliced array with only x node coordinates
#y=jdat1[:,1] # sliced array with only y node coordinates
for i
in range(0,nom):
print ("Member: ", i+1)
x1 = jdat1[mdat1[i,0]-1,0]
y1 = jdat1[mdat1[i,0]-1,1]
x2 = jdat1[mdat1[i,1]-1,0]
y2 = jdat1[mdat1[i,1]-1,1]
print
("joint",mdat1[i,0],":","x1:",x1,"y1:",y1)
print
("joint",mdat1[i,1],":","x2:",x2,"y2:",y2)
# Element length
L = np.sqrt((x2- x1)**2 + (y2-y1)**2)
mprop[i,0] = L
print("L:",L,"mm")
# Cosine and sine of element angle
l = (x2 - x1) / L
m = (y2 - y1) / L
mprop[i,1] = l
mprop[i,2] = m
print("l:",l)
print("m:",m)
#AEL
AEL = mdat4[i] * mdat3[i] / L
mprop[i,3] = AEL
print("AE/L:",AEL,"kN mm^2")
# Local stiffness matrix (4x4)
K = AEL * np.array([[ l*l , l*m ,-l*l ,-l*m],
[ l*m , m*m ,-l*m ,-m*m],
[-l*l ,-l*m , l*l , l*m],
[-l*m ,-m*m , l*m , m*m]])
print(K)
print("mcode",mcode[i,0],mcode[i,1],mcode[i,2],mcode[i,3])
for j in range(0,4):
for k in range(0,4):
if mcode[i,j] != 0 and mcode[i,k] !=
0:
x = mcode[i,j]
y = mcode[i,k]
print("x,y:",x,y)
print("K[",j+1,k+1,"] =",K[j,k])
Ksys[x-1,y-1] =
Ksys[x-1,y-1] + K[j,k]
print("Ksys : \n",Ksys,"\n")
#
Return global stiffness matrix
return Ksys , mprop
#Step
6---------------------------------------------------------------------------
# Define functions
for the Gobal Displacement
def cholesky(Ksys):
print("Step 6: Compute for Gloabal
Displacement-----------------------\n")
print("Ksys : \n",Ksys)
n =
len(Ksys)
L =
np.zeros((n, n))
for i
in range(n):
for j in range(i+1):
if i == j:
L[i, j] = np.sqrt(Ksys[i, i] -
np.sum(L[i, :]**2))
else:
L[i, j] = (Ksys[i, j] - np.sum(L[i,
:] * L[j, :])) / L[j, j]
print("L:",'\n',L)
return L
def
forward_substitution(L, Qsys):
print("Qsys:",'\n',Qsys)
#Get
number of rows
n =
L.shape[0]
#Allocating space for the solution vector
y =
np.zeros_like(Qsys, dtype=np.double);
#Here
we perform the forward-substitution.
#Initializing with the first row.
y[0]
= Qsys[0] / L[0, 0]
#Looping over rows in reverse (from the bottom up),
#starting with the second to last row, because the
#last
row solve was completed in the last step.
for i
in range(1, n):
y[i] = (Qsys[i] - np.dot(L[i,:i], y[:i])) / L[i,i]
print("y:",'\n',y)
return y
def
back_substitution(U, y):
print("U:",'\n',U)
#Number of rows
n =
U.shape[0]
#Allocating space for the solution vector
Dsys
= np.zeros_like(y, dtype=np.double);
#Here
we perform the back-substitution.
#Initializing with the last row.
Dsys[-1] = y[-1] / U[-1, -1]
#Looping over rows in reverse (from the bottom up),
#starting with the second to last row, because the
#last
row solve was completed in the last step.
for i
in range(n-2, -1, -1):
Dsys[i] = (y[i] - np.dot(U[i,i:], Dsys[i:])) / U[i,i]
print("Dsys:",'\n',Dsys,'\n')
return Dsys
# Step
7---------------------------------------------------------------------------
# Define functions
for Member Forces & Reactions
def
Forces(jdat1,mdat1,mcode,mprop,NDF):
print("Step 7: Compute for Member Forces &
Reactions-----------------------\n")
noj =
jdat1.shape[0]
nom =
mdat1.shape[0]
f =
np.zeros((nom))
R =
np.zeros((noj,2))
D =
np.zeros((4))
for i
in range (0,nom):
print ("Member: ", i+1)
for j in range (0,4):
if mcode[i,j] == 0:
D[j]=0
else:
D[j] = Dsys[mcode[i,j]-1]
print("D",j,":",D[j],"mm")
print("AEL =",mprop[i,3])
f[i] = mprop[i,3]*((D[2]-D[0])*mprop[i,1]+(D[3]-D[1])*mprop[i,2])
print("member forces =",f[i],"kN")
Fx1 = -f[i]*mprop[i,1]
Fy1 = -f[i]*mprop[i,2]
Fx2 = f[i]*mprop[i,1]
Fy2 = f[i]*mprop[i,2]
print("Fx1=",Fx1,"kN")
print("Fy1=",Fy1,"kN")
print("Fy2=",Fx2,"kN")
print("Fy2=",Fy2,"kN\n")
R[mdat1[i,0]-1,0] = R[mdat1[i,0]-1,0] + Fx1
R[mdat1[i,0]-1,1] = R[mdat1[i,0]-1,1] + Fy1
R[mdat1[i,1]-1,0] = R[mdat1[i,1]-1,0] + Fx2
R[mdat1[i,1]-1,1] = R[mdat1[i,1]-1,1] + Fy2
for m
in range(0,noj):
print ("Joint: ", m+1)
for n in range(0,2):
if NDF[m,n]!=0:
R[m,n] = 0
print("Reaction
x",m+1,":",R[m,0],"kN")
print("Reaction
y",m+1,":",R[m,1],"kN")
print ("\n")
return f,R
## RUN
# Call for the
title and disclaimer
title_disclaimer()
# Step 1 - Compute
for correctness of in terms of numbers
check(jdat1,jdat3,jdat4,jdat5,mdat1,
mdat3, mdat4)
# Step 2 - Compute
for NDF and TDOF
NDF,TDOF =
DOF(jdat1,jdat3)
# Step 3 - Compute
for Qsys
Qsys =
Q(jdat1,jdat4,jdat5, NDF,TDOF)
# Step 4 - Compute
for Mcode
mcode, HBW =
MCODE(mdat1)
# Step 5 - Compute
for stiffness matrix
Ksys , mprop =
create_stiffness_matrix(jdat1, mdat1, mdat3, mdat4,TDOF,mcode)
# Step 6 - Compute
for Gloabal Displacement
L = cholesky(Ksys)
y =
forward_substitution(L, Qsys)
U = np.transpose(L)
Dsys =
back_substitution(U, y)
#Step 7: Compute
for Member Forces and Reaction
f,R =
Forces(jdat1,mdat1,mcode,mprop,NDF)