Linear Static Structural Analysis

LinearStaticStructuralAnalysis.py

# ------------------------------------------------------------------------
# The following Python code is implemented by Professor Terje Haukaas at
# the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without any form of warranty.
# Please see the Programming note for how to get started, and notice
# that you must copy certain functions into the file terjesfunctions.py
# ------------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt
import terjesfunctions as fcns


# ------------------------------------------------------------------------
# INPUT
# ------------------------------------------------------------------------

# Constants for European IPE220 steel cross-section [N, mm]
E = 200000.0
A = 3340.0
I = 27700000.0
q = 100.0
meter = 1000.0

# Nodes (x, y)
nodes = np.array([(0.0,     0.0),
                  (2*meter, 0.0),
                  (4*meter, 0.0),
                  (5*meter, 0.0)])

# Frame elements (node1, node2, E, I, A, q, nu, beta)
frameElements = np.array([(1, 2, E, I, A, q),
                          (2, 3, E, I, A, q),
                          (3, 4, E, I, A, q)])

# Boundary conditions (node, 0=free, 1=fixed)
constraints = np.array([(1, 1, 1, 0),
                        (3, 0, 1, 0)])

# ------------------------------------------------------------------------
# LOOP OVER ELEMENTS, ASSEMBLE THE STIFFNESS MATRIX Ka
# ------------------------------------------------------------------------

# Number of DOFs
numDOFsPerNodeInStructure = 3
numNodes = nodes.shape[0]
totalNumDOFs = numDOFsPerNodeInStructure * numNodes

# Count elements
numTrussElements = 0
numFrameElements = 0
numQuad4Elements = 0
if 'trussElements' in locals():
    numTrussElements = trussElements.shape[0]
if 'frameElements' in locals():
    numFrameElements = frameElements.shape[0]
if 'quad4Elements' in locals():
    numQuad4Elements = quad4Elements.shape[0]

# Set axial forces for P-delta analysis
if 'trussElements' in locals():
    Ntruss = np.zeros(numTrussElements)
if 'frameElements' in locals():
    Nframe = np.zeros(numFrameElements)


# Initialize the big stiffness matrix
Ka = np.zeros((totalNumDOFs, totalNumDOFs))

# Assemble truss elements
KlStorageTruss = []
TlgStorageTruss = []
for i in range(numTrussElements):

    # Information about element
    E = trussElements[i, 2]
    A = trussElements[i, 3]

    # Information about nodes
    node1 = int(trussElements[i, 0])
    node2 = int(trussElements[i, 1])
    nodeList = np.array([node1, node2])
    nodalcoords = np.array([nodes[node1-1, 0], nodes[node1-1, 1],
                            nodes[node2-1, 0], nodes[node2-1, 1]])

    # Element length and orientation
    delta_x = nodalcoords[2] - nodalcoords[0]
    delta_y = nodalcoords[3] - nodalcoords[1]
    L = np.sqrt(delta_x ** 2 + delta_y ** 2)

    # Get Kl and Tlg
    [Kl, Kgeometric, Tlg] = fcns.getTrussStiffness2D(E, A, L, delta_x, delta_y, Ntruss[i])
    Kl_total = Kl - Kgeometric
    KlStorageTruss.append(Kl_total)
    TlgStorageTruss.append(Tlg)

    # Calculate Kg = Tlg^T * Kl * Tlg
    Kg = (np.transpose(Tlg).dot(Kl_total)).dot(Tlg)

    # Assemble
    numDOFsPerNodeInElement = 2
    Ka = fcns.assembleStiffnessMatrix(Ka, Kg, nodeList, numDOFsPerNodeInElement, numDOFsPerNodeInStructure)

# Assemble frame elements
KlStorageFrame = []
TlgStorageFrame = []
for i in range(numFrameElements):

    # Information about element
    E = frameElements[i, 2]
    I = frameElements[i, 3]
    A = frameElements[i, 4]

    # Information about nodes
    node1 = int(frameElements[i, 0])
    node2 = int(frameElements[i, 1])
    nodeList = np.array([node1, node2])
    nodalcoords = np.array([nodes[node1-1, 0], nodes[node1-1, 1], nodes[node2-1, 0], nodes[node2-1, 1]])

    # Element length and orientation
    delta_x = nodalcoords[2] - nodalcoords[0]
    delta_y = nodalcoords[3] - nodalcoords[1]
    L = np.sqrt(delta_x ** 2 + delta_y ** 2)

    # Get Kl and Tlg
    if len(frameElements[i]) > 7:
        nu = frameElements[i, 6]
        beta = frameElements[i, 7]
    else:
        nu = 0.0
        beta = 0.0
    [Kl, Kgeometric, Tlg] = fcns.getFrameStiffness2D(E, nu, I, A, beta, L, delta_x, delta_y, Nframe[i])
    Kl_total = Kl - Kgeometric
    KlStorageFrame.append(Kl_total)
    TlgStorageFrame.append(Tlg)

    # Calculate Kg = Tlg^T * Kl * Tlg
    Kg = (np.transpose(Tlg).dot(Kl_total)).dot(Tlg)

    # Assemble
    numDOFsPerNodeInElement = 3
    Ka = fcns.assembleStiffnessMatrix(Ka, Kg, nodeList, numDOFsPerNodeInElement, numDOFsPerNodeInStructure)


# Assemble Quad4 elements
for i in range(numQuad4Elements):

    # Information about element
    E = quad4Elements[i, 4]
    nu = quad4Elements[i, 5]
    t = quad4Elements[i, 6]
    ps = quad4Elements[i, 7]

    # Information about nodes
    node1 = int(quad4Elements[i, 0])
    node2 = int(quad4Elements[i, 1])
    node3 = int(quad4Elements[i, 2])
    node4 = int(quad4Elements[i, 3])
    nodeList = np.array([node1, node2, node3, node4])
    nodalcoords = np.array([nodes[node1-1, 0], nodes[node1-1, 1], nodes[node2-1, 0], nodes[node2-1, 1],
                            nodes[node3-1, 0], nodes[node3-1, 1], nodes[node4-1, 0], nodes[node4-1, 1]])

    x = np.array([nodalcoords[0], nodalcoords[2], nodalcoords[4], nodalcoords[6]])
    y = np.array([nodalcoords[1], nodalcoords[3], nodalcoords[5], nodalcoords[7]])

    Kg = fcns.getQuad4Stiffness2D(E, nu, t, ps, x, y)

    # Assemble
    numDOFsPerNodeInElement = 2
    Ka = fcns.assembleStiffnessMatrix(Ka, Kg, nodeList, numDOFsPerNodeInElement, numDOFsPerNodeInStructure)


# ------------------------------------------------------------------------
# ASSEMBLE THE LOAD VECTOR Fa
# ------------------------------------------------------------------------

# Below we set some rough scaling factors for the BMD and SFD for frame members
if 'frameElements' in locals():
    Mmax = 0.0
    Vmax = 0.0
    node1 = int(frameElements[0, 0])
    node2 = int(frameElements[0, 1])
    x1 = nodes[node1 - 1, 0]
    y1 = nodes[node1 - 1, 1]
    x2 = nodes[node2 - 1, 0]
    y2 = nodes[node2 - 1, 1]
    lengthOfFirstElement = np.sqrt((x2-x1)**2 + (y2-y1)**2)

# Create the big load vector
Fa = np.zeros(totalNumDOFs)

# Add nodal loads
if 'loads' in locals():

    numLoads = loads.shape[0]
    for i in range(numLoads):
        loadValue = loads[i, 1] * loads[i, 3]
        Fa[int((loads[i, 0]-1)*numDOFsPerNodeInStructure+(loads[i,  2]-1))] = loadValue

        # Update scaling factors for BMD and SFD
        if 'frameElements' in locals():
            if abs(loadValue) > Vmax:
                Vmax = abs(loadValue)
            if abs(loadValue*lengthOfFirstElement) > Mmax:
                Mmax = abs(loadValue*lengthOfFirstElement)

# Add distributed element loads
if 'frameElements' in locals():
    for i in range(numFrameElements):

        # Get the element load
        if len(frameElements[i]) > 5:
            q = frameElements[i, 5]
        else:
            q = 0.0

        # Proceed only if there is distributed load on this element
        if q != 0.0:

            # Calculate element length
            node1 = int(frameElements[i, 0])
            node2 = int(frameElements[i, 1])
            nodeList = np.array([node1, node2])
            nodalcoords = np.array([nodes[node1-1, 0], nodes[node1-1, 1], nodes[node2-1, 0], nodes[node2-1, 1]])
            delta_x = nodalcoords[2] - nodalcoords[0]
            delta_y = nodalcoords[3] - nodalcoords[1]
            L = np.sqrt(delta_x ** 2 + delta_y ** 2)

            # Local force vector
            Flocal = np.zeros(6)
            Flocal[1] = -q*L/2.0
            Flocal[2] = q*L**2/12.0
            Flocal[4] = -q*L/2.0
            Flocal[5] = -q*L**2/12.0

            # Update scaling factors for BMD and SFD
            if abs(q*L/2.0) > Vmax:
                Vmax = abs(q*L/2.0)
            if abs(q*L**2/12.0) > Mmax:
                Mmax = abs(q*L**2/12.0)

            # Transform it into the global configuration
            Fglobal = (np.transpose(TlgStorageFrame[i])).dot(Flocal)

            # Place it into the force vector for the structure
            numNodesInElement = len(nodeList)
            for j in range(numNodesInElement):
                for m in range(numDOFsPerNodeInElement):
                    Fa[(nodeList[j] - 1) * numDOFsPerNodeInStructure + m] += Fglobal[j * numDOFsPerNodeInElement + m]

# ------------------------------------------------------------------------
# INTRODUCE BOUNDARY CONDITIONS & LOCK ZERO-STIFFNESS DOFs
# ------------------------------------------------------------------------

# DOFs locked by input
dofStatus = np.zeros(totalNumDOFs)
nfix = constraints.shape[0]
for k in range(nfix):
    fixedNode = int(constraints[k, 0])
    for m in range(numDOFsPerNodeInStructure):
        if constraints[k, m + 1] == 1:
            dofStatus[(fixedNode - 1) * numDOFsPerNodeInStructure + m] = 1

# DOFs without stiffness
for i in range(totalNumDOFs):
    if Ka[i, i] == 0:
        dofStatus[i] = 1

numFreeDOFs = int(totalNumDOFs - np.sum(dofStatus))

Kf = np.zeros((numFreeDOFs, numFreeDOFs))
Ff = np.zeros(numFreeDOFs)
iCounter = 0
jCounter = 0
for i in range(totalNumDOFs):

    if dofStatus[i] == 0:

        jCounter = 0

        for j in range(i + 1):

            if dofStatus[j] == 0:
                Kf[iCounter, jCounter] = Ka[i, j]
                Kf[jCounter, iCounter] = Ka[i, j]
                jCounter += 1

        Ff[iCounter] = Fa[i]
        iCounter += 1


# ------------------------------------------------------------------------
# SOLVE THE SYSTEM OF EQUILIBRIUM EQUATIONS AND GET ua
# ------------------------------------------------------------------------
uf = np.linalg.solve(Kf, Ff)

# Get ua by simply adding zero numbers to the array
ua = np.zeros(totalNumDOFs)
iCounter = 0
for i in range(totalNumDOFs):
    if dofStatus[i] == 0:
        ua[i] = uf[iCounter]
        iCounter += 1

# Ouput
print('\n'"Maximum nodal displacement: %.6fmm" % np.max(np.absolute(ua)))


# ------------------------------------------------------------------------
# RECOVER LOCAL ELEMENT FORCES
# ------------------------------------------------------------------------

# Recover truss element forces
for i in range(numTrussElements):

    # From ua to ug
    node1 = int(trussElements[i, 0])
    node2 = int(trussElements[i, 1])
    nodeList = np.array([node1, node2])
    numNodesInElement = len(nodeList)
    numDOFsPerNodeInElement = 2
    ug = np.zeros(numNodesInElement*numDOFsPerNodeInElement)
    for j in range(numNodesInElement):
        for m in range(numDOFsPerNodeInElement):
            ug[j*numDOFsPerNodeInElement+m] = ua[(nodeList[j]-1)*numDOFsPerNodeInStructure+m]

    # From ug to ul
    ul = TlgStorageTruss[i].dot(ug)

    # From ul to Fl
    Fl = KlStorageTruss[i].dot(ul)

    # Get axial force for P-delta analysi (compression positive)
    Ntruss[i] = Fl[0]


# Recover frame element forces
for i in range(numFrameElements):

    # From ua to ug
    node1 = int(frameElements[i, 0])
    node2 = int(frameElements[i, 1])
    nodeList = np.array([node1, node2])
    numNodesInElement = len(nodeList)
    numDOFsPerNodeInElement = 3
    ug = np.zeros(numNodesInElement * numDOFsPerNodeInElement)

    for j in range(numNodesInElement):
        for m in range(numDOFsPerNodeInElement):
            ug[j * numDOFsPerNodeInElement + m] = ua[(nodeList[j] - 1) * numDOFsPerNodeInStructure + m]

    # From ug to ul
    ul = TlgStorageFrame[i].dot(ug)

    # From ul to Fl
    Fl = KlStorageFrame[i].dot(ul)

    # Get axial force for P-delta analysis (compression positive)
    Nframe[i] = Fl[0]


# ------------------------------------------------------------------------
# PRINT AXIAL FORCES IN PREPARATION FOR P-DELTA ANALYSIS
# ------------------------------------------------------------------------

# Notice that compression is positive in this context
if 'trussElements' in locals():
    print('\n'"Axial force in truss members: [", end=" ")
    for i in range(numTrussElements):
        if i == numTrussElements-1:
            print(Ntruss[i], "]")
        else:
            print(Ntruss[i], ",", end=" ")

if 'frameElements' in locals():
    print('\n'"Axial force in frame members: [", end=" ")
    for i in range(numFrameElements):
        if i == numFrameElements-1:
            print(Nframe[i], "]")
        else:
            print(Nframe[i], ",", end=" ")


# ------------------------------------------------------------------------
# PLOT DEFORMED STRUCTURE
# ------------------------------------------------------------------------

# Find outer plot limits
xmin = 0
xmax = 0
ymin = 0
ymax = 0

for i in range(numNodes):

    if nodes[i, 0] < xmin:
        xmin = nodes[i, 0]

    elif nodes[i, 0] > xmax:
        xmax = nodes[i, 0]

    elif nodes[i, 1] < ymin:
        ymin = nodes[i, 1]

    elif nodes[i, 1] > ymax:
        ymax = nodes[i, 1]

xmin = xmin - (xmax - xmin) * 0.2
ymin = ymin - (ymax - ymin) * 0.2
xmax = xmax + (xmax - xmin) * 0.2
ymax = ymax + (ymax - ymin) * 0.2

# Find scaling factor for displacement
x_disp_max = 0.0
y_disp_max = 0.0

for i in range(numNodes):
    dofx = i * numDOFsPerNodeInStructure + 0
    dofy = i * numDOFsPerNodeInStructure + 1

    if abs(ua[dofx]) > x_disp_max:
        x_disp_max = abs(ua[dofx])

    if abs(ua[dofy]) > y_disp_max:
        y_disp_max = abs(ua[dofy])

charact_dim = np.sqrt((xmax-xmin)**2 + (ymax-ymin)**2)

if x_disp_max == 0.0 and y_disp_max == 0.0:
    scale_factor = 0.05 * charact_dim / (0.001 * charact_dim)
    print('\n'"Note: Scaling of displaced shape is probably off because no node displaced")
elif x_disp_max > y_disp_max:
    scale_factor = 0.05 * charact_dim / x_disp_max
else:
    scale_factor = 0.05 * charact_dim / y_disp_max

# Also set scaling factors for the BMD and SFD
if 'frameElements' in locals():
    BMDscaling = 0.05 * charact_dim / Mmax
    SFDscaling = 0.05 * charact_dim / Vmax

# Add max displacement margin around plot
xmin = xmin - x_disp_max * scale_factor
ymin = ymin - y_disp_max * scale_factor
xmax = xmax + x_disp_max * scale_factor
ymax = ymax + y_disp_max * scale_factor

# Check if the structure is 1D
if (ymax - ymin) == 0:
    ymin_new = ymin - y_disp_max * scale_factor * 1.2
    ymax = ymin + y_disp_max * scale_factor * 1.2
    ymin = ymin_new

elif (xmax - xmin) == 0:
    xmin_new = xmin - x_disp_max * scale_factor * 1.2
    xmax = xmin + x_disp_max * scale_factor * 1.2
    xmin = xmin_new

# Make the plotting area square
xrange = xmax-xmin
yrange = ymax-ymin

if xrange > yrange:
    ymin -= 0.5*(xrange-yrange)
    ymax += 0.5*(xrange-yrange)
else:
    xmin -= 0.5*(yrange-xrange)
    xmax += 0.5*(yrange-xrange)

xcoord = np.zeros(2)
ycoord = np.zeros(2)

xdisp = np.zeros(2)
ydisp = np.zeros(2)

xmoment = np.zeros(2)
ymoment = np.zeros(2)

xshear = np.zeros(2)
yshear = np.zeros(2)

if 'frameElements' in locals():
    if len(frameElements[0]) > 7:
        print('\n'"Note: BMD and SFD are NOT plotted when shear deformation is included")

for i in range(numTrussElements):
    node1 = int(trussElements[i, 0])
    node2 = int(trussElements[i, 1])
    xcoord[0] = nodes[node1-1, 0] + ua[(node1 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
    ycoord[0] = nodes[node1-1, 1] + ua[(node1 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor
    xcoord[1] = nodes[node2-1, 0] + ua[(node2 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
    ycoord[1] = nodes[node2-1, 1] + ua[(node2 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor

    plt.plot(xcoord, ycoord, 'ks-')

for i in range(numFrameElements):

    # Get end nodes for this element
    node1 = int(frameElements[i, 0])
    node2 = int(frameElements[i, 1])

    # Get end-point coordinates for this element
    x1 = nodes[node1-1, 0]
    y1 = nodes[node1-1, 1]
    x2 = nodes[node2-1, 0]
    y2 = nodes[node2-1, 1]

    # Calculate element length and direction cosines
    delta_x = x2-x1
    delta_y = y2-y1
    L = np.sqrt(delta_x**2 + delta_y**2)
    cx = delta_x/L
    cy = delta_y/L

    # Get EI for this element
    E = frameElements[i, 2]
    I = frameElements[i, 3]
    EI = E*I

    # Get the element load
    if len(frameElements[i]) > 5:
        q = -frameElements[i, 5]
    else:
        q = 0.0

    # Get local displacements
    nodeList = np.array([node1, node2])
    numNodesInElement = len(nodeList)
    numDOFsPerNodeInElement = 3
    ug = np.zeros(numNodesInElement * numDOFsPerNodeInElement)
    for j in range(numNodesInElement):
        for m in range(numDOFsPerNodeInElement):
            ug[j * numDOFsPerNodeInElement + m] = ua[(nodeList[j] - 1) * numDOFsPerNodeInStructure + m]
    ul = TlgStorageFrame[i].dot(ug)

    # Rename the displacements for clarity
    u1 = ul[0]
    u2 = ul[1]
    u3 = ul[2]
    u4 = ul[3]
    u5 = ul[4]
    u6 = ul[5]

    # Coefficients in the solution to the differential equation
    C1 = -q * L / (12.0 * EI) - 2.0 * u5 / L ** 3 - u3 / L ** 2 + 2.0 * u2 / L ** 3 - u6 / L ** 2
    C2 = q * L * L / (24.0 * EI) + 3.0 * u5 / L ** 2 + 2.0 * u3 / L - 3.0 * u2 / L ** 2 + u6 / L
    C3 = -u3
    C4 = u2

    # Divide the element into "discretization" number of parts
    discretization = 10
    xi = 0.0
    for j in range(discretization):

        # First point of this segment of the element
        x = xi*L
        xA = x1 + (x2 - x1) * xi
        yA = y1 + (y2 - y1) * xi

        # Calculate displacements and bending moment at this point
        u = u1 + xi * (u4 - u1)
        w = q * x**4 / (24.0 * EI) + C1 * x**3 + C2 * x**2 + C3 * x + C4
        M = -(q * x ** 2 / 2.0 + EI * 6.0 * C1 * x + 2.0 * EI * C2)
        V = q * x + EI * 6.0 * C1

        # Rotate and scale displacements, BMD, and SFD into the global coordinate system: global = R * local
        xAdisp = xA + (cx * u - cy * w) * scale_factor
        yAdisp = yA + (cy * u + cx * w) * scale_factor
        xA_BMD = xA - cy * M * BMDscaling
        yA_BMD = yA + cx * M * BMDscaling
        xA_SFD = xA - cy * V * SFDscaling
        yA_SFD = yA + cx * V * SFDscaling

        # Increment xi to get to the second point
        xi += 1.0 / (float(discretization))

        # Second point of this segment of the element
        x = xi*L
        xB = x1 + (x2 - x1) * xi
        yB = y1 + (y2 - y1) * xi

        # Displacement and bending moment at this point
        u = u1 + xi * (u4 - u1)
        w = q * x**4 / (24.0 * EI) + C1 * x**3 + C2 * x**2 + C3 * x + C4
        M = -(q * x ** 2 / 2.0 + EI * 6.0 * C1 * x + 2.0 * EI * C2)
        V = q * x + EI * 6.0 * C1

        # Rotate and scale displacements, BMD, and SFD into the global coordinate system: global = R * local
        xBdisp = xB + (cx * u - cy * w) * scale_factor
        yBdisp = yB + (cy * u + cx * w) * scale_factor
        xB_BMD = xB - cy * M * BMDscaling
        yB_BMD = yB + cx * M * BMDscaling
        xB_SFD = xB - cy * V * SFDscaling
        yB_SFD = yB + cx * V * SFDscaling

        # Plotting for this segment
        xcoord[0] = xA
        xcoord[1] = xB
        ycoord[0] = yA
        ycoord[1] = yB
        plt.plot(xcoord, ycoord, color='silver')

        xdisp[0] = xAdisp
        xdisp[1] = xBdisp
        ydisp[0] = yAdisp
        ydisp[1] = yBdisp
        plt.plot(xdisp, ydisp, 'k-')

        # Only plot BMD and SFD if there isn't shear deformation, because the differential
        # equation solution used above does not include that effect
        if len(frameElements[i]) < 7:

            if j==0:
                xmoment[0] = xA
                xmoment[1] = xA_BMD
                ymoment[0] = yA
                ymoment[1] = yA_BMD
                plt.plot(xmoment, ymoment, 'r-')

                xshear[0] = xA
                xshear[1] = xA_SFD
                yshear[0] = yA
                yshear[1] = yA_SFD
                plt.plot(xshear, yshear, 'b-')

            xmoment[0] = xA_BMD
            xmoment[1] = xB_BMD
            ymoment[0] = yA_BMD
            ymoment[1] = yB_BMD
            plt.plot(xmoment, ymoment, 'r-')

            xshear[0] = xA_SFD
            xshear[1] = xB_SFD
            yshear[0] = yA_SFD
            yshear[1] = yB_SFD
            plt.plot(xshear, yshear, 'b-')

            if j==discretization-1:
                xmoment[0] = xB_BMD
                xmoment[1] = xB
                ymoment[0] = yB_BMD
                ymoment[1] = yB
                plt.plot(xmoment, ymoment, 'r-')

                xshear[0] = xB_SFD
                xshear[1] = xB
                yshear[0] = yB_SFD
                yshear[1] = yB
                plt.plot(xshear, yshear, 'b-')

for i in range(numQuad4Elements):

    node1 = int(quad4Elements[i, 0])
    node2 = int(quad4Elements[i, 1])
    node3 = int(quad4Elements[i, 2])
    node4 = int(quad4Elements[i, 3])

    xcoord[0] = nodes[node1-1, 0] + ua[(node1 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
    ycoord[0] = nodes[node1-1, 1] + ua[(node1 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor
    xcoord[1] = nodes[node2-1, 0] + ua[(node2 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
    ycoord[1] = nodes[node2-1, 1] + ua[(node2 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor
    plt.plot(xcoord, ycoord, 'ks-')

    xcoord[0] = nodes[node3-1, 0] + ua[(node3 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
    ycoord[0] = nodes[node3-1, 1] + ua[(node3 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor
    xcoord[1] = nodes[node2-1, 0] + ua[(node2 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
    ycoord[1] = nodes[node2-1, 1] + ua[(node2 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor
    plt.plot(xcoord, ycoord, 'ks-')

    xcoord[0] = nodes[node3-1, 0] + ua[(node3 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
    ycoord[0] = nodes[node3-1, 1] + ua[(node3 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor
    xcoord[1] = nodes[node4-1, 0] + ua[(node4 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
    ycoord[1] = nodes[node4-1, 1] + ua[(node4 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor
    plt.plot(xcoord, ycoord, 'ks-')

    xcoord[0] = nodes[node1-1, 0] + ua[(node1 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
    ycoord[0] = nodes[node1-1, 1] + ua[(node1 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor
    xcoord[1] = nodes[node4-1, 0] + ua[(node4 - 1) * numDOFsPerNodeInStructure + 0] * scale_factor
    ycoord[1] = nodes[node4-1, 1] + ua[(node4 - 1) * numDOFsPerNodeInStructure + 1] * scale_factor
    plt.plot(xcoord, ycoord, 'ks-')

print('\n'"Colour coding in plot: Grey = Undeformed / Black = Deformed / Red = BMD / Blue = SFD")
plt.axis([xmin, xmax, ymin, ymax])
plt.title("Maximum nodal displacement: %.2fmm" % np.max(np.absolute(ua)))
plt.show()