Dice Arithmetic

Note: It should be run as one page

Choose an operation:

Addition
Subtraction

# of dice:

Attempt: 0

Score: 0

Answer :

Arithmetic

Note: It should be run as one page

Choose an operation:

Addition
Subtraction
Multiplication
Division

Number of digit:

1

Attempt: 0

Score: 0

:

Answer:

Number Game

Note: It should be run as one page

Attempt: 0

Guess a number from 1 to 100:

Mental Calculation of Multiplication by Two-Digit Numbers

n x r = a x c x 100 + (a x d + b x c) x 10 + b x d

let
n = a x 10 + b ; First Two-Digit Number
r = c x 10 + d ; Second Two-Digit Number

example

44 x 51 = (4 x 5 x 100 = 2000) + ((4 x 1 + 4 x 5) x 10 = 240) + (4 x 1 = 4) = 2244

52 x 34 = (5 x 3 x 100 = 1500) + ((5 x 4 + 2 x 3) x 10 = 260) + (2 x 4 = 8) = 1768

61 x 72 = (6 x 7 x 100 = 4200) + ((6 x 2 + 1 x 7) x 10 = 190) + (1 x 2 = 2) = 4392

91 x 24 = (9 x 2 x 100 = 1800) + ((9 x 4 + 1 x 2) x 10 = 380) + (1 x 4 = 4) = 2184

Mental Calculation of n^2 Base 100

n^2 = (n+(n-100))x100+(n-100)^2

From

n^2 = (n-100)^2 + 200n - 100^2 <<< (n^2 - 200n + 100^2) + 200n - 100^2

example

105^2 = ((105 + 5) x 100 = 11000) + ( 5 x 5 = 25) = 11025

107^2 = ((107 + 7) x 100 = 11400) + ( 7 x 7 = 49) = 11449

111^2 = ((111 + 11) x 100 = 12200) + (11 x 11 = 121) = 12321

97^2 = (( 97 - 3) x 100 = 9400) + (-3 x -3 = 9) = 9409

e Learning

AppSheet e Learning

Overview:

An interactive app for math, vocabulary, memorization enthusiasts.

Challenges users with two-digit math problems, Word games, & Flashcard.

Includes addition, subtraction, Multiply, Dictionary, and User Define Flash Card Data.

Random questions.

Scoring system.

Immediate feedback.

Customization options.

How to Use:

Launch the app.

Start the quiz.

Solve math problems, word problem, and memorization.

Check your results.

Aim for a high score!

Open AppSheet Link

Math Quiz 2 Digits

AppSheet Math Quiz 2 Digits

Overview:

An interactive app for math enthusiasts.  Challenges users with two-digit math problems.

Includes addition, subtraction, and Multiply. 

Random questions.  Scoring system. Immediate feedback. Time System.

How to Use:

Launch the app. Start the quiz. Solve math problems. Check your results.

Aim for a high score!

Open AppSheet Link

Propped Using Double Integration Method with Mathcad

by Ernel Solis Filipinas (ernel_filipinas@yahoo.com).

 

https://drive.google.com/file/d/1PtNurT-ZXgc4HmASNSCz4cDi2q9de8eP/view?usp=drive_link

 

(For educational reference only)

Simply Supported Using Double Integration Method with Mathcad

by Ernel Solis Filipinas (ernel_filipinas@yahoo.com).

 

https://drive.google.com/file/d/1IzM_4rCl6s6R4vTUtwJZ_xGuX2hLnJcR/view?usp=drive_link

 

(For educational reference only)

Cantilever Using Double Integration Method with Mathcad

by Ernel Solis Filipinas (ernel_filipinas@yahoo.com).

 

https://drive.google.com/file/d/1qrjNU0HFYj6NookR9zsmCAfW4q3j07Hk/view?usp=drive_link

 

(For educational reference only)

2D BEAM ANALYSIS (DIRECT STIFFNESS METHOD) PYTHON PROGRAMMING

2D BEAM ANALYSIS v4 first created as FRAME ANALYSIS on 24 July 2009 using Excel VBA and updated on 10 July 2024 using Phyton by Ernel Solis Filipinas (ernel_filipinas@yahoo.com).

 

https://colab.research.google.com/drive/1wuFzXgiJxAPUaEyDAd3oQ2KnrpKOHKtZ?usp=sharing

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 m

jdat1 = np.array([

    [0,0],

    [10,0]

    ])

# Define node/joint type [] where 1-Hinge, 2-Roller x, 3-Roller y, 4-Free, 5-Fixed Joints

jdat2 = np.array([[5],[2]])

# Define  Moment Joint Load  in kNm

# External Force (+ve counter counter clockwise, -ve clockwise)

# Fixed End force (inverse of External Force)

jdat5 = np.array([[-83.33333],[83.33333]])

# Define Vertical Joint Load  in kN

# External Force (+ve going up, -ve going down)

# Fixed End force (inverse of External Force)

jdat6 = np.array([[-50],[-50]])

# Define element/member incidences (joint end1, joint end2)

mdat1 = np.array([[1,2]])

# Define Modulus of Elasticity (kPa)

mdat2 = np.array([[1]])

# Define Moment of Inertia (m4)

mdat3 = np.array([[1]])

# Define Fixed End Moment for beam kNm # Make sure it is as per member incidence

#(+ve counter counter clockwise, -ve clockwise)

mdat4 = np.array([[83.33333,-83.33333]])

# Define Fixed End Shear for beam kNm

#(+ve going up, -ve going down)

mdat5 = np.array([[50,50]])

#FEM Link

#https://engineering.purdue.edu/~ce474/Docs/Fixed%20End%20Moments.pdf

# Scale

dscale = 5000     #Deformation Scale

ascale = 200      #Arrow Scale

#----------------------------------------------------------------------------------

#----------------------------------------------------------------------------------

## CALCULATION FUNCTION

# Step 0---------------------------------------------------------------------------

def title_disclaimer():

    print("2D BEAM ANALYSIS v4 first created as FRAME ANALYSIS in 09 January 2021 using Excel VBA")

    print("10 July 2024 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,jdat2,jdat5,jdat6,mdat1, mdat2, mdat3, mdat4, mdat5):

    print("Step 1: Checking number of joint and number of member-----------------------\n")

    noj = jdat1.shape[0]

    nojdat5 = jdat5.shape[0]

    nojdat6 = jdat6.shape[0]

    print ("Number of joints =" ,noj)

    if noj == nojdat5 and noj == nojdat6:

        print("The number joints and joint datum is OK")

    else:

        sys.exit("Number of joints not equal for jdat data")

    nom = mdat1.shape[0]

    nonmdat2 = mdat2.shape[0]

    nonmdat3 = mdat3.shape[0]

    nonmdat4 = mdat4.shape[0]

    nonmdat5 = mdat5.shape[0]

    print ("Number of members " ,nom)

    if nom == nonmdat2 and nom==nonmdat3 and nom==nonmdat4 and nom==nonmdat5:

        print("The number members and member datum is OK\n")

    else:

        sys.exit("Number of members not equal to nomdat data")

#Step 2---------------------------------------------------------------------------

# Define functions for NDF, TDOF

def DOF(jdat1,jdat2):

    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 jdat2[i] == 1: #Hinge

            NDF[i,0] = 0

            NDF[i,1] = TDOF + 1

            TDOF+=1

        elif jdat2[i] == 2: #Roller x

            NDF[i,0] = 0

            NDF[i,1] = TDOF + 1

            TDOF+=1

        elif jdat2[i] == 3: #Roller y

            NDF[i,0] = TDOF + 1

            NDF[i,1] = TDOF + 2

            TDOF+=2

        elif jdat2[i] == 4: #Free

            NDF[i,0] = TDOF + 1

            NDF[i,1] = TDOF + 2

            TDOF+=2

        elif jdat2[i] == 5: #Fixed

            NDF[i,0] = 0

            NDF[i,1] = 0

        else:

            print("Please check the joint data")

        print("NDF Array Y = ", NDF[i,0],"NDF Array M = ", NDF[i,1])

    print ("TDOF = ",TDOF,"\n")

    return NDF, TDOF

#Step 3---------------------------------------------------------------------------

# Define functions for Mcode

def MCODE(mdat1,NDF, TDOF):

    print("Step 4: Compute for Mcode-----------------------\n")

    nom = mdat1.shape[0]

    mcode = np.zeros(shape=(nom,4),dtype = int) # for beam only

    MDIFF = 0

    for i in range(0,nom):

        print("Member:",i+1)

        print("Member Incidence",mdat1[i,0],mdat1[i,1])

        #mcode[i,] = 0 # For beam only

        mcode[i,0] = NDF[mdat1[i,0]-1,0]

        mcode[i,1] = NDF[mdat1[i,0]-1,1]

        #mcode[i,] = 0 # For beam only

        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 4---------------------------------------------------------------------------

# Define functions for Qsys

def Q(jdat1,jdat5,jdat6, NDF,TDOF,MCODE):

    print("Step 3: Compute for Qsys-----------------------\n")

    noj = jdat1.shape[0]

    Qsys = np.zeros(shape=(TDOF),dtype = float)

    nom = mdat1.shape[0]

    F = np.zeros(shape=(noj),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] = jdat6[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],"kNm")

    return Qsys

#Step 5---------------------------------------------------------------------------

# Define functions for Ksys

def create_stiffness_matrix(jdat1, mdat1, mdat2, mdat3,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,6))

    # 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,"m")

        # 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)

        #EI

        EIL = mdat2[i] * mdat3[i]/L

        mprop[i,3] = EIL.item()

        mprop[i,4] = (6*EIL.item())/(L**1)

        mprop[i,5] = (12*EIL.item())/(L**2)

        print("EI/L:",EIL,"kNm/rad")

        k11 = 12/(L**2)

        k12 = 6/(L**1)

        # Local stiffness matrix (4x4)

        K = EIL * np.array([[ k11 , k12 ,-k11 , k12],

                           [ k12 , 4 ,-k12 , 2],

                           [-k11 ,-k12 , k11 ,-k12],

                           [ k12 , 2 ,-k12 , 4]

                           ])

        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]

    # Allocate space for the solution vector

    Dsys = np.zeros_like(y, dtype=np.double)

    # Initialize with the last row

    Dsys[-1] = y[-1] / U[-1, -1]

    # Loop over rows in reverse order

    for i in range(n - 2, -1, -1):

        Dsys[i] = (y[i] - np.dot(U[i, i + 1:], Dsys[i + 1:])) / U[i, i]

    print("Dsys:",'\n',Dsys,'\n')

    noj = jdat1.shape[0]

    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 z",m+1,":",Dsys[NDF[m,n]-1],"m")

                if n==1:

                    print("Rotation y",m+1,":",Dsys[NDF[m,n]-1],"rad\n")

            if NDF[m,n]==0:

                if n==0:

                    print("Displacement z",m+1,":",0,"m")

                if n==1:

                    print("roation y",m+1,":",0,"rad\n")

    return Dsys

# Step 7---------------------------------------------------------------------------

# Define functions for Member Forces & Reactions

def Forces(jdat1,mdat1,mcode,mprop,NDF,Dsys):

    print("Step 7: Compute for Member Forces & Reactions-----------------------\n")

    noj = jdat1.shape[0]

    nom = mdat1.shape[0]

    MDOF = np.zeros((nom,2))

    FEM = np.zeros((nom,2))

    M = np.zeros((nom,2))

    RDOF = np.zeros((nom,2))

    FER = np.zeros((nom,2))

    R = np.zeros((nom,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("D1",D[0],"m")

        print("D2",D[1],"rad")

        print("D3",D[2],"m")

        print("D4",D[3],"rad")

        print("EI/L =",mprop[i,3])

        MDOF[i,0] = mprop[i,3]*(4*D[1]+2*D[3])+mprop[i,4]*(D[0]-D[2])

        MDOF[i,1] = mprop[i,3]*(2*D[1]+4*D[3])+mprop[i,4]*(D[0]-D[2])

        print("member Moment force end 1 due to DOF = ", MDOF[i,0],"kNm")

        print("member Moment force end 2 due to DOF = ", MDOF[i,1],"kNm")

        FEM[i,0] = mdat4[i,0]

        FEM[i,1] = mdat4[i,1]

        print("member Moment force end 1 due to FEM = ", FEM[i,0],"kNm")

        print("member Moment force end 2 due to FEM = ", FEM[i,1],"kNm")

        M[i,0] = MDOF[i,0]+FEM[i,0]

        M[i,1] = MDOF[i,1]+FEM[i,1]

        print("member Moment force end 1 Final = ", M[i,0],"kNm")

        print("member Moment force end 2 Final = ", M[i,1],"kNm\n")

        RDOF[i,0] = mprop[i,4]*(D[1]+D[3])+mprop[i,5]*(D[0]-D[2])

        RDOF[i,1] = -mprop[i,4]*(D[1]+D[3])+mprop[i,5]*(D[2]-D[0])

        print("member Vertical force end 1 due to DOF = ", RDOF[i,0],"kN")

        print("member Vertical force end 2 due to DOF = ", RDOF[i,1],"kN")

        FER[i,0] = mdat5[i,0]

        FER[i,1] = mdat5[i,1]

        print("member Vertica force end 1 due to FEM = ", FER[i,0],"km")

        print("member Vertica force end 2 due to FEM = ", FER[i,1],"km")

        R[i,0] = RDOF[i,0]+FER[i,0]

        R[i,1] = RDOF[i,1]+FER[i,1]

        print("member Vertica force end 1 Final = ", R[i,0],"kN")

        print("member Vertica force end 2 Final = ", R[i,1],"kN\n")

    return M, R

# Step 8---------------------------------------------------------------------------

# Define functions for Graph

def TGraph(jdat1,jdat2,mdat1,mcode):

    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 = 'none', 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

        plt.plot(xp,yp, color='green', linestyle='solid', linewidth = 2, marker = ',', ms = 4)

    #Adding support

    for m in range (0,noj):

        Sx1 = jdat1[m,0]

        Sx2 = jdat1[m,1]

        if jdat2[m] == 1:

            plt.plot(Sx1,Sx2, color='Blue', marker = '^', ms = 10)

        elif jdat2[m] == 2:

            plt.plot(Sx1,Sx2, color='Blue', marker = 'o', ms = 10)

        elif jdat2[m] == 3:

            plt.plot(Sx1,Sx2, color='Blue', marker = 'o', ms = 10)

        elif jdat2[m] == 4:

            plt.plot(Sx1,Sx2, color='Blue', marker = ',', ms = 10)

        elif jdat2[m] == 5:

            plt.plot(Sx1,Sx2, color='Blue', marker = 's', ms = 10)

        plt.annotate(m+1, (jdat1[m,0]+ascale/10, jdat1[m,1]+ascale/10))

    # 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 Beam 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,jdat2,jdat5,jdat6,mdat1, mdat2, mdat3,mdat4, mdat5)

# Step 2 - Compute for NDF and TDOF

NDF,TDOF = DOF(jdat1,jdat2)

# Step 3 - Compute for Mcode

mcode, HBW = MCODE(mdat1,NDF, TDOF)

# Step 4 - Compute for Qsys

Qsys = Q(jdat1,jdat5,jdat6, NDF,TDOF,MCODE)

# Step 5 - Compute for stiffness matrix

Ksys , mprop = create_stiffness_matrix(jdat1, mdat1, mdat2, mdat3,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

M,R = Forces(jdat1,mdat1,mcode,mprop,NDF,Dsys)

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/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)