How to do pushover analysis using Python

print("==========================")

from openseespy.opensees import *


import os, sys
import numpy as np
import math as math
import matplotlib.pyplot as plt


print("Starting RCFrameGravity example")

wipe()
# Create ModelBuilder (with two-dimensions and 3 DOF/node)
model('basic', '-ndm', 2, '-ndf', 3)

# Create nodes
# ------------

# Set parameters for overall model geometry
height = 1300.0

# Create nodes
#    tag, X, Y
node(1, 0.0, 0.0)
node(2, 0.0, 0.0)
node(3, 0.0, height)

# Fix supports at base of columns
#   tag, DX, DY, RZ
fix(1, 1, 1, 1)
fix(2, 1, 1, 0)

#%%
# Define materials
uniaxialMaterial('Hysteretic', 1, 0.006, 75514606.25, 0.03, 75514606.25*1.1, 
                 -0.006, -75514606.25, -0.03, -75514606.25*1.1, 1.0, 1.0, 0.0, 0.0, 0.7)
#%%
# Define elements
# ----------------------

# Geometry of column elements
#                tag

geomTransf('PDelta', 1)

# Number of integration points along length of element
np = 10

# Lobatto integratoin
beamIntegration('Lobatto', 1, 1, np)

# Create the coulumns 
element('zeroLength',  1,   1,   2,  '-mat',1, '-dir',6)
bc = 250.0
dc = 250.0
A  = bc*dc
Iy = bc*dc**3/12
Ec = 26692.696
element('elasticBeamColumn',  2,   2,   3,  A,  Ec,  Iy, 1)

# import vfo rendering module
import vfo.vfo as vfo

# render the model after defining all the nodes and elements
# vfo.plot_model(show_eletags="yes")

# plot mode shape
# vfo.plot_modeshape(modenumber=3)

#%%

# Define gravity loads
# --------------------

#  a parameter for the axial load
P = 200000.0

# Create a Plain load pattern with a Linear TimeSeries
timeSeries('Linear', 1)
pattern('Plain', 1, 1)

# Create nodal loads at nodes 3 
#    nd  FX,  FY, MZ
# load(3, -P, 0.0, 0.0)
load(3, 0.0, -P, 0.0)

# ------------------------------
# End of model generation
# ------------------------------


# ------------------------------
# Start of analysis generation
# ------------------------------

# Create the system of equation, a sparse solver with partial pivoting
system('ProfileSPD')

# Create the constraint handler, the transformation method
constraints('Transformation')

# Create the DOF numberer, the reverse Cuthill-McKee algorithm
numberer('RCM')

# Create the convergence test, the norm of the residual with a tolerance of
# 1e-12 and a max number of iterations of 50
test('NormDispIncr', 1.0e-12, 50, 3)

# Create the solution algorithm, a Newton-Raphson algorithm
algorithm('Newton')

# Create the integration scheme, the LoadControl scheme using steps of 0.1
integrator('LoadControl', 0.1)

# Create the analysis object
analysis('Static')

# ------------------------------
# End of analysis generation
# ------------------------------


# ------------------------------
# Finally perform the analysis
# ------------------------------

# perform the gravity load analysis, requires 50 steps to reach the load level
analyze(50)

# Print out the state of nodes 3 
# print node 3 

# Print out the state of element 1 and 2
# print ele 1 2

u3 = nodeDisp(3, 1)

results = open('results.out', 'a+')

results.close()

print("==========================")

#%%

# ----------------------------------------------------
# Start of Model Generation & Initial Gravity Analysis
# ----------------------------------------------------

# Set the gravity loads to be constant & reset the time in the domain
loadConst('-time', 0.0)

# ----------------------------------------------------
# End of Model Generation & Initial Gravity Analysis
# ----------------------------------------------------


# ----------------------------------------------------
# Start of additional modelling for lateral loads
# ----------------------------------------------------

# Define lateral loads
# --------------------

# Set some parameters
H = 200000.0  # Reference lateral load

# Set lateral load pattern with a Linear TimeSeries
pattern('Plain', 2, 1)

# Create nodal loads at nodes 3
#    nd    FX  FY  MZ
load(3, H, 0.0, 0.0)

# ----------------------------------------------------
# End of additional modelling for lateral loads
# ----------------------------------------------------

# ----------------------------------------------------
# Start of modifications to analysis for push over
# ----------------------------------------------------

# Set some parameters
dU = 1.0  # Displacement increment

# Change the integration scheme to be displacement control
#                             node dof init Jd min max
integrator('DisplacementControl', 3, 1, dU, 1, dU, dU)

# ----------------------------------------------------
# End of modifications to analysis for push over
# ----------------------------------------------------

# ------------------------------
# Start of recorder generation
# ------------------------------

# Stop the old recorders by destroying them
# remove recorders

# Create a recorder to monitor nodal displacements
# recorder Node -file node3.out -time -node 3 -dof 1 2 3 disp

# Create a recorder to monitor element forces in columns
# recorder EnvelopeElement -file ele12.out -time -ele 1 2 forces

# --------------------------------
# End of recorder generation
# ---------------------------------

# ------------------------------
# Finally perform the analysis
# ------------------------------

# Set some parameters
maxU = 250.0  # Max displacement
currentDisp = 0.0
ok = 0

test('NormDispIncr', 1.0e-12, 1000)
algorithm('ModifiedNewton', '-initial')


topdispl = [] #empty list to record top node displacement
ndispl = {1: [], 2: [], 3: []}
eforce = {(1,1): [], (1,2):[], (1,3):[], (1,4):[], (1,5):[], (1,6):[], 
          (2,1): [], (2,2):[], (2,3):[], (2,4):[], (2,5):[], (2,6):[]}

while ok == 0 and currentDisp < maxU:

    ok = analyze(1)
    print('currentDisp =',currentDisp)

    # if the analysis fails try initial tangent iteration
    if ok != 0:
        print("modified newton failed")
        break
    # print "regular newton failed .. lets try an initail stiffness for this step"
    # test('NormDispIncr', 1.0e-12,  1000)
    # # algorithm('ModifiedNewton', '-initial')
    # ok = analyze(1)
    # if ok == 0:
    #     print "that worked .. back to regular newton"

    # test('NormDispIncr', 1.0e-12,  10)
    # algorithm('Newton')

    currentDisp = nodeDisp(3, 1)
    
    #record each step
    topdispl.append(nodeDisp(3,1))
    for ntag in [1,2,3]:
        ndispl[ntag].append(nodeDisp(ntag,1))
        
    for etag in [1,2]:
        eforce[etag,1].append(eleForce(etag,1))
        eforce[etag,2].append(eleForce(etag,2))
        eforce[etag,3].append(eleForce(etag,3))
        eforce[etag,4].append(eleForce(etag,4))
        eforce[etag,5].append(eleForce(etag,5))
        eforce[etag,6].append(eleForce(etag,6))
       

xpoints=ndispl[3]
ypoints=eforce[2,1]
plt.plot(xpoints, ypoints)

    
results = open('results.out', 'a+')

results.close()

# Print the state at node 3
# print node 3

print("==========================")