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/10ekgFQ7Ct8fJmWWrXz27NTp_7wCggbN8
For educational purpose only
----------------------------------------------------------------------------------------------------------------------------
import math
import numpy as np
import matplotlib.pyplot as plt
import sys
## INPUT DATA
# Define joint coordinate [x,y] in mm
jdat1 = np.array([
[0,0],
[1000,0],
[0,1000],
[1000,1000]
])
# Define node/joint type [] where 1-Hinge, 2-Roller x, 3-Roller y, 4-Free Joint
jdat3 = np.array([[1],[2],[4],[4]])
# Define external horizontal load in kN (+ve going right, -ve going left)
jdat4 = np.array([[100],[0],[1000],[0]])
# Define external vertical load in kN (+ve going up, -ve going down)
jdat5 = np.array([[100],[0],[0],[0]])
# Define element/member incidences (joint end1, joint end2)
mdat1 = np.array([[1,2],[3,4],[1,3],[2,4],[2,3],[1,4]])
# Define Area in mm2
mdat3 = np.array([[2500],[2500],[2500],[2500],[2500],[2500]])
# Define Modulus of Elasticity (GPa)
mdat4 = np.array([[200],[200],[200],[200],[200],[200]])
# Scale
dscale = 500 #Deformation Scale
ascale = 200 #Arrow Scale
#----------------------------------------------------------------------------------
#----------------------------------------------------------------------------------
## CALCULATION FUNCTION
# Step 0---------------------------------------------------------------------------
def title_disclaimer():
print("2D TRUSS ANALYSIS v4.3 first created on 14 June 2001 using QBasic, ")
print("then updated on 09 December 2023 using Phyton")
print("and provided with basic graph on 07 July 2024 ")
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 organization 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,number of member, joint data, & member data\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:
sys.exit("Number of joints not equal to number of vertical loads")
else:
sys.exit("Number of joints not equal to number of horizontal loads ")
else:
sys.exit("Number of joints not equal to number of joint type")
nom = mdat1.shape[0]
nonmdat3 = mdat3.shape[0]
nonmdat4 = mdat4.shape[0]
print ("")
print ("Number of members " ,nom)
if nom == nonmdat3:
if nom == nonmdat4:
print("The number members and member datum is OK\n")
else:
sys.exit("Number of members not equal to number of member areas\n")
else:
sys.exit("Number of members not equal to number of member modulus\n")
print("\nJoint Input Data\n")
for i in range (0,noj):
print("Joint: ", i+1)
print("x coordanate : ",jdat1[i,0])
print("y coordanate : ",jdat1[i,1])
print("Joint Type : ",jdat3[i])
print("Horizontal Load : ",jdat4[i])
print("Vertical Load : ",jdat5[i],"\n")
print("\nMemeber Input Data\n")
for i in range (0,nom):
print("Member:",i+1)
print("Member Incidence",mdat1[i,0],mdat1[i,1])
print("Member Area : ",mdat3[i])
print("Member Modulus : ",mdat4[i],"\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 = float)
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].item()
print("Qsys",NDF[i,0],"=", Qsys[NDF[i,0]-1],"kN")
if j==1:
Qsys[NDF[i,1]-1] = jdat5[i].item()
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.item()
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]
noj = jdat1.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')
for m in range(0,noj):
print ("Joint: ", m+1)
for n in range(0,2):
if NDF[m,n]!=0:
if n==0:
print("Displacement x",m+1,":",Dsys[NDF[m,n]-1],"mm")
if n==1:
print("Displacement y",m+1,":",Dsys[NDF[m,n]-1],"mm\n")
if NDF[m,n]==0:
if n==0:
print("Displacement x",m+1,":",0,"mm")
if n==1:
print("Displacement y",m+1,":",0,"mm\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))
mdis = np.zeros((nom))
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]
mdis = ((D[2]-D[0])*mprop[i,1]+(D[3]-D[1])*mprop[i,2])
print ("member displacement=",mdis,"mm")
print("Fx1=",Fx1,"kN")
print("Fy1=",Fy1,"kN")
print("Fx2=",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
Rx = np.round(R[m,0]-jdat4[m], decimals=6)
Ry = np.round(R[m,1]-jdat5[m], decimals=6)
print("Reaction x",m+1,":",Rx,"kN")
print("Reaction y",m+1,":",Ry,"kN\n")
return f,R
# Step 8---------------------------------------------------------------------------
# Define functions for Graph
def TGraph(jdat1,jdat3,mdat1,NDF,Dsys,Qsys):
noj = jdat1.shape[0]
nom = mdat1.shape[0]
dist = np.zeros((4))
Qsys_x = np.zeros(shape=(TDOF),dtype = float)
Qsys_y = np.zeros(shape=(TDOF),dtype = float)
for i in range(0,nom):
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]
# x axis values
xp =np.array([x1,x2])
# corresponding y axis values
yp =np.array([y1,y2])
# plotting the points
plt.plot(xp,yp, color='green', linestyle='solid', linewidth = 2, marker = 'o', ms = 4)
for j in range (0,4):
if mcode[i,j] == 0:
dist[j]=0
else:
dist[j] = Dsys[mcode[i,j]-1]
dx1 = x1 + dist[0]*dscale
dy1 = y1 + dist[1]*dscale
dx2 = x2 + dist[2]*dscale
dy2 = y2 + dist[3]*dscale
# x axis values for displacement
dxp =np.array([dx1,dx2])
# corresponding y axis values jfor displacement
dyp =np.array([dy1,dy2])
# plotting the points of displacement
plt.plot(dxp,dyp, color='Cyan', linestyle='dashed', linewidth = 1, marker = 'o', ms = 4)
#Adding support
for m in range (0,noj):
if jdat3[m] == 1:
Sx1 = jdat1[m,0]
Sx2 = jdat1[m,1]
plt.plot(Sx1,Sx2, color='Blue', marker = '^', ms = 10)
elif jdat3[m] == 2:
Sx1 = jdat1[m,0]
Sx2 = jdat1[m,1]
plt.plot(Sx1,Sx2, color='Blue', marker = 'o', ms = 10)
elif jdat3[m] == 3:
Sx1 = jdat1[m,0]
Sx2 = jdat1[m,1]
plt.plot(Sx1,Sx2, color='Blue', marker = 'o', ms = 10)
plt.annotate(m+1, (jdat1[m,0]+ascale/10, jdat1[m,1]+ascale/10))
for n in range (0,2):
if n==0:
if jdat4[m] > 0:
plt.arrow(x=jdat1[m,0]+ascale/5, y= jdat1[m,1], dx=ascale, dy=0, width=ascale/50,head_width = ascale/10,facecolor='red',edgecolor='red')
plt.annotate(jdat4[m],xy=(jdat1[m,0]+ascale/5+ascale,jdat1[m,1]+ascale/10),horizontalalignment='left',color = 'red')
elif jdat4[m] < 0:
plt.arrow(x=jdat1[m,0]-ascale/5, y= jdat1[m,1], dx=-ascale, dy=0, width=ascale/50,head_width = ascale/10,facecolor='red',edgecolor='red')
plt.annotate(jdat4[m],xy=(jdat1[m,0]-ascale/5-ascale,jdat1[m,1]+ascale/10),horizontalalignment='left',color = 'red')
if n==0:
if jdat5[m] > 0:
plt.arrow(x=jdat1[m,0], y= jdat1[m,1]+ascale/5, dx=0, dy=ascale, width=ascale/50,head_width = ascale/10,facecolor='red',edgecolor='red')
plt.annotate(jdat5[m],xy=(jdat1[m,0]+ascale/10,jdat1[m,1]+ascale/5+ascale),horizontalalignment='left',color = 'red')
elif jdat5[m] < 0:
plt.arrow(x=jdat1[m,0], y= jdat1[m,1]+ascale/5+ascale, dx=0, dy=-ascale, width=ascale/50,head_width = ascale/10,facecolor='red',edgecolor='red')
plt.annotate(jdat5[m],xy=(jdat1[m,0]+ascale/10,jdat1[m,1]+ascale/5+ascale),horizontalalignment='left',color = 'red')
# naming the x axis
plt.xlabel('x - axis')
# naming the y axis
plt.ylabel('y - axis')
# giving a title to my graph
plt.title('2D Truss Matrix Stiffness Analysis with Graph \n by Ernel S. Filipinas in 2024')
# function to show the plot
plt.show()
## 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)