terjehaukaas/mdofBucklingAnalysis.py

## mdofBucklingAnalysis.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 [N, m, s, kg]
# ------------------------------------------------------------------------

# Constants
E = 2.0e+11
I = 8.09e-5
A = 0.0091
storeyHeight = 4.0
bayWidth = 6.0

# Nodes (x, y)
nodes = np.array([(0,           0),
                  (bayWidth,    0),
                  (2*bayWidth,  0),
                  (0,           storeyHeight),
                  (bayWidth,    storeyHeight),
                  (2*bayWidth,  storeyHeight),
                  (0,           2*storeyHeight),
                  (bayWidth,    2*storeyHeight),
                  (2*bayWidth,  2*storeyHeight)])

# Frame elements (node1, node2, E, I, A)
frameElements = np.array([(1, 4, E, I, A),
                          (2, 5, E, I, A),
                          (3, 6, E, I, A),
                          (4, 5, E, I, A),
                          (5, 6, E, I, A),
                          (4, 7, E, I, A),
                          (5, 8, E, I, A),
                          (6, 9, E, I, A),
                          (7, 8, E, I, A),
                          (8, 9, E, I, A)])

# Nodal load multipliers for the buckling load (node, value, dof)
bucklingLoadMultiplier = np.array([(4, -1, 2),
                                   (5, -1, 2),
                                   (6, -1, 2),
                                   (7, -1, 2),
                                   (8, -1, 2),
                                   (9, -1, 2)])

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


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

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

# Count elements
numFrameElements = frameElements.shape[0]

# Set axial force in elements
Nframe = np.zeros(numFrameElements)

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

# 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
    nu = 0
    beta = 0
    [Kl_elastic, Kl_geometric, Tlg] = fcns.getFrameStiffness2D(E, nu, I, A, beta, L, delta_x, delta_y, Nframe[i])
    KlStorageFrame.append(Kl_elastic)
    TlgStorageFrame.append(Tlg)

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

    # Assemble
    numDOFsPerNodeInElement = 3
    Ka_elastic = fcns.assembleStiffnessMatrix(Ka_elastic, Kg_elastic, nodeList, numDOFsPerNodeInElement, numDOFsPerNodeInStructure)


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

# Notice that only point loads are included for now
Fa = np.zeros(totalNumDOFs)

numBucklingLoadMultipliers = bucklingLoadMultiplier.shape[0]

for i in range(numBucklingLoadMultipliers):
    Fa[(bucklingLoadMultiplier[i, 0]-1)*numDOFsPerNodeInStructure+(bucklingLoadMultiplier[i,  2]-1)] \
        = bucklingLoadMultiplier[i, 1]


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

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

Kf_elastic = 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_elastic[iCounter, jCounter] = Ka_elastic[i, j]
                Kf_elastic[jCounter, iCounter] = Ka_elastic[i, j]
                jCounter += 1

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

# ------------------------------------------------------------------------
# SOLVE THE SYSTEM OF EQUILIBRIUM EQUATIONS AND GET ua
# ------------------------------------------------------------------------
uf = np.linalg.solve(Kf_elastic, 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


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

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


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

# Initialize the big stiffness matrices
Ka_geometric = np.zeros((totalNumDOFs, totalNumDOFs))

# Assemble frame elements
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
    nu = 0
    beta = 0
    [Kl_elastic, Kl_geometric, Tlg] = fcns.getFrameStiffness2D(E, nu, I, A, beta, L, delta_x, delta_y, Nframe[i])


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

    # Assemble
    numDOFsPerNodeInElement = 3
    Ka_geometric = fcns.assembleStiffnessMatrix(Ka_geometric, Kg_geometric, nodeList, numDOFsPerNodeInElement, numDOFsPerNodeInStructure)

# ------------------------------------------------------------------------
# INTRODUCE BOUNDARY CONDITIONS
# ------------------------------------------------------------------------

Kf_geometric = np.zeros((numFreeDOFs, 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_geometric[iCounter, jCounter] = Ka_geometric[i, j]
                Kf_geometric[jCounter, iCounter] = Ka_geometric[i, j]

                jCounter += 1

        iCounter += 1


# ------------------------------------------------------------------------
# EIGENVALUE ANALYSIS
# ------------------------------------------------------------------------

# 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

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

# Here we need scipy, to solve the *generalized* eigenvalue problem [K-omega^2*M]{u} = 0
import scipy

# Determine eigenvalues
allEigenvalues, allEigenvectors = scipy.linalg.eig(Kf_elastic, Kf_geometric)

eigenvalues = []
eigenvectors = []
for k in range(len(allEigenvalues)):

    if allEigenvalues[k].real == np.inf:

        pass

    else:

        eigenvalues.append(allEigenvalues[k].real)

        modeShape = []

        # Extract the mode
        for j in range(numFreeDOFs):

            modeShape.append(allEigenvectors[j, k])

        eigenvectors.append(modeShape)

# Sort the modes
idx = np.array(eigenvalues).argsort()[::]

# Print summary of buckling loads
print('\n'"Found the following", len(eigenvalues), "buckling loads:")


for i in range(len(eigenvalues)):

    index = idx[i]

    print("P = %.1f kN" % (0.001*eigenvalues[index]))

    # Clear plot
    plt.cla()

    # Get ua by adding zero numbers to the array
    ua = np.zeros(totalNumDOFs)
    mode = eigenvectors[index]
    jCounter = 0
    for j in range(totalNumDOFs):
        if dofStatus[j] == 0:
            ua[j] = mode[jCounter]
            jCounter += 1

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

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

    scale_factor = 1/np.max(np.abs(mode))

    for j in range(numFrameElements):

        # Get end nodes for this element
        node1 = int(frameElements[j, 0])
        node2 = int(frameElements[j, 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[j, 2]
        I = frameElements[j, 3]
        EI = E * I

        # Get local displacements
        nodeList = np.array([node1, node2])
        numNodesInElement = len(nodeList)
        numDOFsPerNodeInElement = 3
        ug = np.zeros(numNodesInElement * numDOFsPerNodeInElement)
        for n in range(numNodesInElement):
            for m in range(numDOFsPerNodeInElement):
                ug[n * numDOFsPerNodeInElement + m] = ua[(nodeList[n] - 1) * numDOFsPerNodeInStructure + m]
        ul = TlgStorageFrame[j].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 = -2.0 * u5 / L ** 3 - u3 / L ** 2 + 2.0 * u2 / L ** 3 - u6 / L ** 2
        C2 = 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 = C1 * x ** 3 + C2 * x ** 2 + C3 * x + C4

            # 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

            # 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 = ul[0] + xi * (ul[3] - ul[0])
            w = C1 * x ** 3 + C2 * x ** 2 + C3 * x + C4

            # 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

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

    plt.ion()
    plt.axis([xmin, xmax, ymin, ymax])
    plt.title("Buckling Mode for $P_{cr}=%.1f kN$" % (0.001*eigenvalues[index]))
    plt.show()
    plt.pause(2)
	# ------------------------------------------------------------------------
	# 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 [N, m, s, kg]
	# ------------------------------------------------------------------------

	# Constants
	E = 2.0e+11
	I = 8.09e-5
	A = 0.0091
	storeyHeight = 4.0
	bayWidth = 6.0

	# Nodes (x, y)
	nodes = np.array([(0, 0),
	(bayWidth, 0),
	(2*bayWidth, 0),
	(0, storeyHeight),
	(bayWidth, storeyHeight),
	(2*bayWidth, storeyHeight),
	(0, 2*storeyHeight),
	(bayWidth, 2*storeyHeight),
	(2bayWidth, 2storeyHeight)])

	# Frame elements (node1, node2, E, I, A)
	frameElements = np.array([(1, 4, E, I, A),
	(2, 5, E, I, A),
	(3, 6, E, I, A),
	(4, 5, E, I, A),
	(5, 6, E, I, A),
	(4, 7, E, I, A),
	(5, 8, E, I, A),
	(6, 9, E, I, A),
	(7, 8, E, I, A),
	(8, 9, E, I, A)])

	# Nodal load multipliers for the buckling load (node, value, dof)
	bucklingLoadMultiplier = np.array([(4, -1, 2),
	(5, -1, 2),
	(6, -1, 2),
	(7, -1, 2),
	(8, -1, 2),
	(9, -1, 2)])

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


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

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

	# Count elements
	numFrameElements = frameElements.shape[0]

	# Set axial force in elements
	Nframe = np.zeros(numFrameElements)

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

	# 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
	nu = 0
	beta = 0
	[Kl_elastic, Kl_geometric, Tlg] = fcns.getFrameStiffness2D(E, nu, I, A, beta, L, delta_x, delta_y, Nframe[i])
	KlStorageFrame.append(Kl_elastic)
	TlgStorageFrame.append(Tlg)

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

	# Assemble
	numDOFsPerNodeInElement = 3
	Ka_elastic = fcns.assembleStiffnessMatrix(Ka_elastic, Kg_elastic, nodeList, numDOFsPerNodeInElement, numDOFsPerNodeInStructure)


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

	# Notice that only point loads are included for now
	Fa = np.zeros(totalNumDOFs)

	numBucklingLoadMultipliers = bucklingLoadMultiplier.shape[0]

	for i in range(numBucklingLoadMultipliers):
	Fa[(bucklingLoadMultiplier[i, 0]-1)*numDOFsPerNodeInStructure+(bucklingLoadMultiplier[i, 2]-1)] \
	= bucklingLoadMultiplier[i, 1]


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

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

	Kf_elastic = 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_elastic[iCounter, jCounter] = Ka_elastic[i, j]
	Kf_elastic[jCounter, iCounter] = Ka_elastic[i, j]
	jCounter += 1

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

	# ------------------------------------------------------------------------
	# SOLVE THE SYSTEM OF EQUILIBRIUM EQUATIONS AND GET ua
	# ------------------------------------------------------------------------
	uf = np.linalg.solve(Kf_elastic, 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


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

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


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

	# Initialize the big stiffness matrices
	Ka_geometric = np.zeros((totalNumDOFs, totalNumDOFs))

	# Assemble frame elements
	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
	nu = 0
	beta = 0
	[Kl_elastic, Kl_geometric, Tlg] = fcns.getFrameStiffness2D(E, nu, I, A, beta, L, delta_x, delta_y, Nframe[i])


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

	# Assemble
	numDOFsPerNodeInElement = 3
	Ka_geometric = fcns.assembleStiffnessMatrix(Ka_geometric, Kg_geometric, nodeList, numDOFsPerNodeInElement, numDOFsPerNodeInStructure)

	# ------------------------------------------------------------------------
	# INTRODUCE BOUNDARY CONDITIONS
	# ------------------------------------------------------------------------

	Kf_geometric = np.zeros((numFreeDOFs, 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_geometric[iCounter, jCounter] = Ka_geometric[i, j]
	Kf_geometric[jCounter, iCounter] = Ka_geometric[i, j]

	jCounter += 1

	iCounter += 1


	# ------------------------------------------------------------------------
	# EIGENVALUE ANALYSIS
	# ------------------------------------------------------------------------

	# 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

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

	# Here we need scipy, to solve the generalized eigenvalue problem [K-omega^2*M]{u} = 0
	import scipy

	# Determine eigenvalues
	allEigenvalues, allEigenvectors = scipy.linalg.eig(Kf_elastic, Kf_geometric)

	eigenvalues = []
	eigenvectors = []
	for k in range(len(allEigenvalues)):

	if allEigenvalues[k].real == np.inf:

	pass

	else:

	eigenvalues.append(allEigenvalues[k].real)

	modeShape = []

	# Extract the mode
	for j in range(numFreeDOFs):

	modeShape.append(allEigenvectors[j, k])

	eigenvectors.append(modeShape)

	# Sort the modes
	idx = np.array(eigenvalues).argsort()[::]

	# Print summary of buckling loads
	print('\n'"Found the following", len(eigenvalues), "buckling loads:")


	for i in range(len(eigenvalues)):

	index = idx[i]

	print("P = %.1f kN" % (0.001*eigenvalues[index]))

	# Clear plot
	plt.cla()

	# Get ua by adding zero numbers to the array
	ua = np.zeros(totalNumDOFs)
	mode = eigenvectors[index]
	jCounter = 0
	for j in range(totalNumDOFs):
	if dofStatus[j] == 0:
	ua[j] = mode[jCounter]
	jCounter += 1

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

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

	scale_factor = 1/np.max(np.abs(mode))

	for j in range(numFrameElements):

	# Get end nodes for this element
	node1 = int(frameElements[j, 0])
	node2 = int(frameElements[j, 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[j, 2]
	I = frameElements[j, 3]
	EI = E * I

	# Get local displacements
	nodeList = np.array([node1, node2])
	numNodesInElement = len(nodeList)
	numDOFsPerNodeInElement = 3
	ug = np.zeros(numNodesInElement * numDOFsPerNodeInElement)
	for n in range(numNodesInElement):
	for m in range(numDOFsPerNodeInElement):
	ug[n * numDOFsPerNodeInElement + m] = ua[(nodeList[n] - 1) * numDOFsPerNodeInStructure + m]
	ul = TlgStorageFrame[j].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 = -2.0 * u5 / L 3 - u3 / L 2 + 2.0 * u2 / L 3 - u6 / L 2
	C2 = 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 = C1 * x ** 3 + C2 * x ** 2 + C3 * x + C4

	# 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

	# 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 = ul[0] + xi * (ul[3] - ul[0])
	w = C1 * x ** 3 + C2 * x ** 2 + C3 * x + C4

	# 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

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

	plt.ion()
	plt.axis([xmin, xmax, ymin, ymax])
	plt.title("Buckling Mode for $P_{cr}=%.1f kN$" % (0.001*eigenvalues[index]))
	plt.show()
	plt.pause(2)