Sunday, December 10, 2023

2D TRUSS ANALYSIS (DIRECT STIFFNESS METHOD) PYTHON PROGRAMMING


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)