PushOversteel.py

import openseespy.opensees as ops
import matplotlib.pyplot as plt
import csv
import os

from Gravity_Analysis import *

print("---------------------------")
print("Gravity Analysis Done.")
print("---------------------------")

refLoad = 1000.0  # N

ops.timeSeries('Linear', 11)
ops.pattern('Plain', 11, 11)

push_direction = 1   # 1, 2 for Pushover in X and Y Direction Respectively
control_node = 5     # node where displacement is read

if push_direction == 2:
    print("Starting Pushover Analysis In Y Direction....")
    control_DOF = 2
elif push_direction == 1:
    print("Starting Pushover Analysis In X Direction....")
    control_DOF = 1
else:
    print("ERROR Pushover Direction.")

if push_direction == 2:
    for tag in master_nodes:
        ops.load(tag, 0.0, refLoad, 0.0, 0.0, 0.0, 0.0)
else:
    for tag in master_nodes:
        ops.load(tag, refLoad, 0.0, 0.0, 0.0, 0.0, 0.0)

dU = 0.05
ops.integrator('DisplacementControl', control_node, control_DOF, dU, 1, dU, dU)

maxDisp = 0.04 * 12600
currentDisp = 0.0
controlNode_disp = []
base_shear = []
ok = 0

ops.test('NormDispIncr', 1.0e-6, 1000)
ops.algorithm('Newton')

# ── Create CSV and write header immediately ──────────────────────────────────
csv_filename = "pushover_results_steel_jacketing.csv"
csv_path = os.path.abspath(csv_filename)
print(f"CSV will be saved to: {csv_path}")

csv_file = open(csv_filename, mode='w', newline='')
writer = csv.writer(csv_file)
writer.writerow(["Drift (%)", "Base Shear (kN)"])   # Write header once
# ─────────────────────────────────────────────────────────────────────────────

temp = 1

try:
    while ok == 0 and currentDisp < maxDisp:

        ok = ops.analyze(1)

        if ok != 0:
            print("Newton failed. Trying different algorithms...")
            algorithms = [
                ('ModifiedNewton', ['-initial']),
                ('NewtonLineSearch', []),
                ('KrylovNewton', []),
                ('BFGS', [])
            ]
            for alg, args in algorithms:
                ops.algorithm(alg, *args)
                if ops.analyze(1) == 0:
                    print(f"Succeeded with {alg}. Back to regular Newton.")
                    ok = 0
                    break
            ops.algorithm('Newton')

        currentDisp = ops.nodeDisp(control_node, control_DOF)
        controlNode_disp.append(currentDisp)

        ops.reactions()
        bShear = 0.0
        for node in floor_1_nodes:
            reaction = -ops.nodeReaction(node, control_DOF)
            bShear += reaction
            ops.nodeResponse(node, control_DOF, 6)
        base_shear.append(bShear / 1.0e3)

        # ── Write each step to CSV immediately ──────────────────────────────
        current_drift = (currentDisp / (40. * ft)) * 100
        writer.writerow([round(current_drift, 6), round(bShear / 1.0e3, 6)])
        csv_file.flush()                 # Force write to disk after every step
        # ────────────────────────────────────────────────────────────────────

        if temp % 20 == 0:
            print(f"Disp : Node {control_node} : {currentDisp:.3f} mm, Base Shear : {(bShear / 1000):.3f} kN")
        temp += 1

finally:
    csv_file.close()     # Always close the file, even if analysis crashes
    print(f"CSV saved: {temp - 1} steps written to '{csv_filename}'")

# ── Summary ──────────────────────────────────────────────────────────────────
if base_shear:
    print(f"Maximum Base Shear = {max(base_shear):.2f} kN")
    max_index = base_shear.index(max(base_shear))
    disp_at_max_base_shear = controlNode_disp[max_index]
    print(f"Displacement at Maximum Base Shear = {disp_at_max_base_shear:.2f} mm")
    drift_at_max = (disp_at_max_base_shear / (40. * ft)) * 100
    print(f"Drift at Maximum Base Shear = {drift_at_max:.2f} %")
else:
    print("WARNING: No data collected. Check if analysis started correctly.")

# ── Plot ─────────────────────────────────────────────────────────────────────
drift = [(d / (40. * ft)) * 100 for d in controlNode_disp]

plt.plot(drift, base_shear, label='Control Node')
plt.title("Pushover Curve")
plt.xlabel("Drift (%)")
plt.xlim(-0.05, 4.05)
plt.ylim(-5, 1.1 * max(base_shear) if base_shear else 100)
plt.ylabel("Base Shear (kN)")
plt.tick_params(direction='in', top=True, right=True)
plt.grid(True, linestyle='--', alpha=0.4)
plt.tight_layout()
plt.show()