## Example of RC column with mixed beam element using OpenSeePy

Forum for asking and answering questions related to use of the OpenSeesPy module

Moderators: silvia, selimgunay, Moderators

Prafullamalla
Posts: 135
Joined: Mon Feb 02, 2015 6:32 pm

### Example of RC column with mixed beam element using OpenSeePy

"""
Created on Thu May 28 16:27:29 2020

@author: praf_malla@hotmail.com, Nepal
@ Corona time, Lockdown
"""
# =============================================================================
# ANALYSIS OF RC COLUMN
# =============================================================================
# Units in N, mm, sec
print("==========================")

from openseespy.postprocessing.Get_Rendering import *
from openseespy.opensees import *
wipe()

import numpy as np
import matplotlib.pyplot as plt
from math import asin, sqrt

print("Starting RC column analysis")
PI = 2 * asin(1.0)

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

# Set parameters for overall model geometry
height = 2000.0 # height of column

# Create nodes
node(1, 0.0, 0.0)
node(2, 0.0, height)

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

# Define cross-section for nonlinear columns
# some parameters
cover = 40 ; # concrete cover
barDia = 20 ; # longitudinal bar dia
As = PI/4*barDia**2 # area of longitudinal bar

# Properties of concrete
# Define parameters of concrete
fc = -30; # Concrete compressive strength (+Tension, -compression)
Ec = 5700*sqrt(-fc); # Modulus of elasticity of concrete
e0 = 0.002; # strain at compressive strength of concrete

# Define parameters of stirrups for confinement
fyh = 415; # yield strength of stirrups
diaStirrup = 8 ; # Diameter of stirrups
(nLY, nLZ) = (2, 3) # no of Legs of stirrup parallel to Y and Z - direction
sh = 150 ; # Spacing of stirrups

# Derived confinement parameters from stirrups
hP = colDepth-2*cover; # Depth of concrete core to oustide of stirrups
bP = colWidth-2*cover; # Width of concrete core to oustide of stirrups
Aso = PI/4*diaStirrup**2; # Area of stirrups
Vstirrup = Aso*(nLY*hP+nLZ*bP)*1000/sh; # Volume of stirrups per 1000 mm
Vcore = hP*bP*1000; # Volume of concrete core in 1000 mm
Rhos = Vstirrup/Vcore; # ratio of volume of hoop steel to the volume of concrete core measured to the outside of hoops.

# Parameters of Confined concrete
Kfc = 1+Rhos*fyh/(-fc); # Ratio of confined to unconfined concrete strength (by Scott)
fc1C = -Kfc*fc; # Compressive strength of confined concrete
eps1C = Kfc*e0; # Strain at confined compressive strength
fc2C = 0.2*fc1C; # stress of confined concrete at ultimate strain

e50uC = (0.0207+eps1C*fc1C)/(fc1C-6.89)
e50hC = 3/4*Rhos*sqrt(hP/sh)
Zc = 0.5/(e50uC+e50hC-eps1C) ; # Descending branch of Confined concrete
eps2C = 0.8*fc1C/Zc+eps1C; # strain at ultimate stress of confined concrete
ftUC = 0.14*fc1C; EtsC = 2*ftUC/eps1C; lambdaC = 0.1

# Parameters of Unconfined concrete
fc1U =-fc ; # Compressive strength of unconfined concrete
eps1U = e0 ; # strain at compressvie strength of unconfined concrete
fc2U = 0.2*fc1U ; # stress of unconfined concrete at ultimate strain

e50uU = (0.0207+eps1U*fc1U)/(fc1U-6.89)
e50hU = 0.0 ; # There is no sirrups i.e., Rhos = 0.0
Zu = 0.5/(e50uU+e50hU-eps1U); # Descending branch of Unconfined concrete
eps2U = 0.8*fc1U/Zu+eps1U ; # Strain at ultimate stress of unconfined concrete
ftUU = 0.14*fc1U; EtsU = 2*ftUU/eps1U; lambdaU = 0.1

# ------------------------------------------
# CONCRETE tag
# Core concrete (confined)
uniaxialMaterial('Concrete02', 1, -fc1C, -eps1C, -fc2C, -eps2C, lambdaC , ftUC, EtsC)
# Cover concrete (unconfined)
uniaxialMaterial('Concrete02', 2, -fc1U, -eps1U, -fc2U, -eps2U, lambdaU , ftUU, EtsU)

print ("Parameters of Confined concrete", -fc1C, -eps1C, -fc2C, -eps2C)
print ("Parameters of UnConfined concrete", -fc1U, -eps1U, -fc2U, -eps2U)

# Plot stress-strain curve of Confined and Unconfined concrete
plt.figure()
sigC = (0, fc1C, fc2C); epsC = (0, eps1C, eps2C)
sigU = (0, fc1U, fc2U); epsU = (0, eps1U, eps2U)
plt.plot(epsC, sigC, linewidth=2.0)
plt.plot(epsU, sigU, linewidth=2.0)
plt.xlabel('strain')
plt.ylabel('compressive stress MPa')

# Longitudinal Reinforcing steel
fy = 400.0; # Yield stress
E = 200000.0; # Young's modulus
Bs = 0.01; # strain hardening ratio
(R0, cR1, cR2) = (18, 0.925, 0.15) # Steel parameters

# Variables derived for formation of section
coverY = colDepth / 2.0
coverZ = colWidth / 2.0
coreY = coverY-cover
coreZ = coverZ-cover
(nfCoreY, nfCoreZ, nfCoverY, nfCoverZ)= (20, 20, 20, 20)
secTag = 1
# tag
uniaxialMaterial('Steel02', 3, fy, E, Bs, R0, cR1, cR2)
section('Fiber', secTag)

# Create the concrete core fibers
patch('quad', 1, nfCoreZ, nfCoreY, -coreY, coreZ, -coreY, -coreZ, coreY, -coreZ, coreY, coreZ,)

# Create the concrete cover patches (right,left,bottom, top)
patch('quad', 2, 2, nfCoverY, -coverY, coverZ, - coreY, coreZ, coreY, coreZ, coverY, coverZ)
patch('quad', 2, 2, nfCoverY, -coreY, -coreZ, -coverY, -coverZ, coverY, -coverZ, coreY, -coreZ)
patch('quad', 2, nfCoverZ, 2, -coverY, coverZ, -coverY, -coverZ, -coreY, -coreZ, -coreY, coreZ)
patch('quad', 2, nfCoverZ, 2, coreY, coreZ, coreY, -coreZ, coverY, -coverZ, coverY, coverZ)

# Create the reinforcing fibers (top, middle, bottom)
layer('straight', 3, 3, As, coreY, coreZ, coreY, -coreZ)
layer('straight', 3, 2, As, 0.0, -coreZ, 0.0, coreZ)
layer('straight', 3, 3, As, -coreY, coreZ, -coreY, -coreZ)

# Define column elements
# ----------------------
geomTag = 1
# Geometry of column elements
# tag
geomTransf('Corotational', geomTag)
np = 4 # Number of integration points along length of element

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

# Create the coulumns using Beam-column elements
# e tag ndI ndJ transfTag integrationTag
#eleType = 'forceBeamColumn'
#element(eleType, 1, 1, 2, 1, 1)

# mixedBeamColumn tag iNode jNode numIntgrPts secTag transftag
eleType = 'mixedBeamColumn'
element(eleType, 1, 1, 2, np, secTag, geomTag, '-integration', 'Lobatto')

# a parameter for the axial load
AxialLoad = 0.1*fc1U*colDepth*colWidth; # 10% of axial capacity of columns

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

# Create nodal loads at nodes 3 & 4
# nd FX, FY, MZ
# ------------------------------
# End of model generation
# ------------------------------

# ------------------------------
# Start of analysis generation
# ------------------------------
# Gravity Analysis
# Create the system of equation, a sparse solver with partial pivoting
Tol = 1.0e-6
system ('BandGeneral')
constraints ('Plain')
numberer ('RCM')
test ('NormDispIncr', Tol, 10, 3)
algorithm ('Newton')
analysis ('Static')
analyze (10)

# Display Model
plot_model()

print("==========================")
print("Gravity Analysis Completed")

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

print("Pushover Analysis Started!")
# Set the gravity loads to be constant & reset the time in the domain
# Set some parameters
# Set lateral load pattern with a Linear TimeSeries
pattern('Plain', 2, 1)

# Create nodal loads at nodes 3 & 4
# nd FX FY MZ

# ----------------------------------------------------
# Start of push over analysis
# ----------------------------------------------------

# Set some parameters
dU = 1 # Displacement increment
dUMin = 1
dUMax = 1

(Tol, MaxIter) = (1.0e-12, 10) # Tolerance and Maximum iteration
(controlNode, controlDOF) =(2, 1) # Control node and control DOF
system ('BandGeneral')
constraints ('Plain')
numberer ('RCM')
test ('NormDispIncr', Tol, MaxIter)
analysis ('Static');

# Change the integration scheme to be displacement control
# node dof init Jd min max
integrator ('DisplacementControl', controlNode, controlDOF, dU, 1, dUMin,dUMax)
algorithm ('Newton')

# Set some parameters
maxU = 140.0 # Max displacement
currentDisp = nodeDisp(controlNode, controlDOF)
ok = 0

# perform the analysis
Nsteps=int(maxU/dU)
import numpy as np
Nsteps=int(maxU/dU) # Number of data to be stored
data=np.zeros( (Nsteps+1, 2) ) # Initial variable to store data
j = 0 # Counter to store data
while (ok == 0 and currentDisp < maxU):
j+=1
analyze(1)
# if the analysis fails try initial tangent iteration
if ok != 0:
print("regular newton failed .. lets try an initail stiffness for this step")
test('NormDispIncr', 1.0e-6, 1000)
algorithm('ModifiedNewton', '-initial')
ok = analyze(1)
algorithm('Newton')

currentDisp = nodeDisp(controlNode, controlDOF)
data[j,0]=currentDisp
plt.plot(data[:,0], data[:,1], linewidth=2.0)
plt.xlabel('Lateral displacement')
plt.show()
wipe

Prafullamalla
Posts: 135
Joined: Mon Feb 02, 2015 6:32 pm

### Re: Example of RC column with mixed beam element using OpenSeePy

Looks, there will be problem on sharing files in this platform as python is indent based program

mhscott
Posts: 444
Joined: Tue Jul 06, 2004 3:38 pm
Location: Corvallis, Oregon USA
Contact:

### Re: Example of RC column with mixed beam element using OpenSeePy

Thanks for sharing the example!

Please try posting again with BBcode.
Do you need extra help with OpenSees? https://courses.silviasbrainery.com/

Michael H. Scott, Ph.D.
OpenSees blog: www.portwooddigital.com

Prafullamalla
Posts: 135
Joined: Mon Feb 02, 2015 6:32 pm

### Re: Example of RC column with mixed beam element using OpenSeePy

Thanks, Prof. Scott. Yeah, it's working. Its been few days, I have started using Python.

"""
Created on Thu May 28 16:27:29 2020

@author: praf_malla@hotmail.com, Nepal
Corona virus time, Lockdown.
Stay home, Keep community safe.
"""
# =============================================================================
# ANALYSIS OF RC COLUMN
# =============================================================================
# Units in N, mm, sec
print("==========================")

from openseespy.postprocessing.Get_Rendering import *
from openseespy.opensees import *
wipe()

import numpy as np
import matplotlib.pyplot as plt
from math import asin, sqrt

print("Starting RC column analysis")
PI = 2 * asin(1.0)

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

# Set parameters for overall model geometry
height = 2000.0 # height of column

# Create nodes
node(1, 0.0, 0.0)
node(2, 0.0, height)

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

# Define cross-section for nonlinear columns
# some parameters
cover = 40 ; # concrete cover
barDia = 20 ; # longitudinal bar dia
As = PI/4*barDia**2 # area of longitudinal bar

# Properties of concrete
# Define parameters of concrete
fc = -30; # Concrete compressive strength (+Tension, -compression)
Ec = 5700*sqrt(-fc); # Modulus of elasticity of concrete
e0 = 0.002; # strain at compressive strength of concrete

# Define parameters of stirrups for confinement
fyh = 415; # yield strength of stirrups
diaStirrup = 8 ; # Diameter of stirrups
(nLY, nLZ) = (2, 3) # no of Legs of stirrup parallel to Y and Z - direction
sh = 150 ; # Spacing of stirrups

# Derived confinement parameters from stirrups
hP = colDepth-2*cover; # Depth of concrete core to oustide of stirrups
bP = colWidth-2*cover; # Width of concrete core to oustide of stirrups
Aso = PI/4*diaStirrup**2; # Area of stirrups
Vstirrup = Aso*(nLY*hP+nLZ*bP)*1000/sh; # Volume of stirrups per 1000 mm
Vcore = hP*bP*1000; # Volume of concrete core in 1000 mm
Rhos = Vstirrup/Vcore; # ratio of volume of hoop steel to the volume of concrete core measured to the outside of hoops.

# Parameters of Confined concrete
Kfc = 1+Rhos*fyh/(-fc); # Ratio of confined to unconfined concrete strength (by Scott)
fc1C = -Kfc*fc; # Compressive strength of confined concrete
eps1C = Kfc*e0; # Strain at confined compressive strength
fc2C = 0.2*fc1C; # stress of confined concrete at ultimate strain

e50uC = (0.0207+eps1C*fc1C)/(fc1C-6.89)
e50hC = 3/4*Rhos*sqrt(hP/sh)
Zc = 0.5/(e50uC+e50hC-eps1C) ; # Descending branch of Confined concrete
eps2C = 0.8*fc1C/Zc+eps1C; # strain at ultimate stress of confined concrete
ftUC = 0.14*fc1C; EtsC = 2*ftUC/eps1C; lambdaC = 0.1

# Parameters of Unconfined concrete
fc1U =-fc ; # Compressive strength of unconfined concrete
eps1U = e0 ; # strain at compressvie strength of unconfined concrete
fc2U = 0.2*fc1U ; # stress of unconfined concrete at ultimate strain

e50uU = (0.0207+eps1U*fc1U)/(fc1U-6.89)
e50hU = 0.0 ; # There is no sirrups i.e., Rhos = 0.0
Zu = 0.5/(e50uU+e50hU-eps1U); # Descending branch of Unconfined concrete
eps2U = 0.8*fc1U/Zu+eps1U ; # Strain at ultimate stress of unconfined concrete
ftUU = 0.14*fc1U; EtsU = 2*ftUU/eps1U; lambdaU = 0.1

# ------------------------------------------
# CONCRETE tag
# Core concrete (confined)
uniaxialMaterial('Concrete02', 1, -fc1C, -eps1C, -fc2C, -eps2C, lambdaC , ftUC, EtsC)
# Cover concrete (unconfined)
uniaxialMaterial('Concrete02', 2, -fc1U, -eps1U, -fc2U, -eps2U, lambdaU , ftUU, EtsU)

print ("Parameters of Confined concrete", -fc1C, -eps1C, -fc2C, -eps2C)
print ("Parameters of UnConfined concrete", -fc1U, -eps1U, -fc2U, -eps2U)

# Plot stress-strain curve of Confined and Unconfined concrete
plt.figure()
sigC = (0, fc1C, fc2C); epsC = (0, eps1C, eps2C)
sigU = (0, fc1U, fc2U); epsU = (0, eps1U, eps2U)
plt.plot(epsC, sigC, linewidth=2.0)
plt.plot(epsU, sigU, linewidth=2.0)
plt.xlabel('strain')
plt.ylabel('compressive stress MPa')

# Longitudinal Reinforcing steel
fy = 400.0; # Yield stress
E = 200000.0; # Young's modulus
Bs = 0.01; # strain hardening ratio
(R0, cR1, cR2) = (18, 0.925, 0.15) # Steel parameters

# Variables derived for formation of section
coverY = colDepth / 2.0
coverZ = colWidth / 2.0
coreY = coverY-cover
coreZ = coverZ-cover
(nfCoreY, nfCoreZ, nfCoverY, nfCoverZ)= (20, 20, 20, 20)
secTag = 1
# tag
uniaxialMaterial('Steel02', 3, fy, E, Bs, R0, cR1, cR2)
section('Fiber', secTag)

# Create the concrete core fibers
patch('quad', 1, nfCoreZ, nfCoreY, -coreY, coreZ, -coreY, -coreZ, coreY, -coreZ, coreY, coreZ,)

# Create the concrete cover patches (right,left,bottom, top)
patch('quad', 2, 2, nfCoverY, -coverY, coverZ, - coreY, coreZ, coreY, coreZ, coverY, coverZ)
patch('quad', 2, 2, nfCoverY, -coreY, -coreZ, -coverY, -coverZ, coverY, -coverZ, coreY, -coreZ)
patch('quad', 2, nfCoverZ, 2, -coverY, coverZ, -coverY, -coverZ, -coreY, -coreZ, -coreY, coreZ)
patch('quad', 2, nfCoverZ, 2, coreY, coreZ, coreY, -coreZ, coverY, -coverZ, coverY, coverZ)

# Create the reinforcing fibers (top, middle, bottom)
layer('straight', 3, 3, As, coreY, coreZ, coreY, -coreZ)
layer('straight', 3, 2, As, 0.0, -coreZ, 0.0, coreZ)
layer('straight', 3, 3, As, -coreY, coreZ, -coreY, -coreZ)

# Define column elements
# ----------------------
geomTag = 1
# Geometry of column elements
# tag
geomTransf('Corotational', geomTag)

np = 4 # Number of integration points along length of element
# Lobatto integration
#beamIntegration('Lobatto', 1, 1, np)

# Create the coulumns using Beam-column elements
# e tag ndI ndJ transfTag integrationTag
#eleType = 'forceBeamColumn'
#element(eleType, 1, 1, 2, 1, 1)

# mixedBeamColumn tag iNode jNode numIntgrPts secTag transftag
eleType = 'mixedBeamColumn'
element(eleType, 1, 1, 2, np, secTag, geomTag, '-integration', 'Lobatto')

# a parameter for the axial load
AxialLoad = 0.1*fc1U*colDepth*colWidth; # 10% of axial capacity of columns

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

# Create nodal loads at nodes 3 & 4
# nd FX, FY, MZ
# ------------------------------
# End of model generation
# ------------------------------

# ------------------------------
# Start of analysis generation
# ------------------------------
# Gravity Analysis
# Create the system of equation, a sparse solver with partial pivoting
Tol = 1.0e-6
system ('BandGeneral')
constraints ('Plain')
numberer ('RCM')
test ('NormDispIncr', Tol, 10, 3)
algorithm ('Newton')
analysis ('Static')
analyze (10)

# Display Model
plot_model()

print("==========================")
print("Gravity Analysis Completed")

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

print("Pushover Analysis Started!")
# Set the gravity loads to be constant & reset the time in the domain
# Set some parameters
# Set lateral load pattern with a Linear TimeSeries
pattern('Plain', 2, 1)

# Create nodal loads at nodes 3 & 4
# nd FX FY MZ

# ----------------------------------------------------
# Start of push over analysis
# ----------------------------------------------------

# Set some parameters
dU = 1 # Displacement increment
dUMin = 1
dUMax = 1

(Tol, MaxIter) = (1.0e-12, 10) # Tolerance and Maximum iteration
(controlNode, controlDOF) =(2, 1) # Control node and control DOF
system ('BandGeneral')
constraints ('Plain')
numberer ('RCM')
test ('NormDispIncr', Tol, MaxIter)
analysis ('Static');

# Change the integration scheme to be displacement control
# node dof init Jd min max
integrator ('DisplacementControl', controlNode, controlDOF, dU, 1, dUMin,dUMax)
algorithm ('Newton')

# Set some parameters
maxU = 140.0 # Max displacement
currentDisp = nodeDisp(controlNode, controlDOF)
ok = 0

# perform the analysis
Nsteps=int(maxU/dU)
import numpy as np
Nsteps=int(maxU/dU) # Number of data to be stored
data=np.zeros( (Nsteps+1, 2) ) # Initial variable to store data
j = 0 # Counter to store data
while (ok == 0 and currentDisp < maxU):
j+=1
analyze(1)
# if the analysis fails try initial tangent iteration
if ok != 0:
print("regular newton failed .. lets try an initail stiffness for this step")
test('NormDispIncr', 1.0e-6, 1000)
algorithm('ModifiedNewton', '-initial')
ok = analyze(1)
algorithm('Newton')

currentDisp = nodeDisp(controlNode, controlDOF)
data[j,0]=currentDisp
plt.plot(data[:,0], data[:,1], linewidth=2.0)
plt.xlabel('Lateral displacement')
plt.show()
wipe

Prafullamalla
Posts: 135
Joined: Mon Feb 02, 2015 6:32 pm

### Re: Example of RC column with mixed beam element using OpenSeePy

Prafullamalla wrote:
Sat May 30, 2020 8:12 am
Thanks, Prof. Scott. Yeah, it's working. Its been few days, I have started using Python.

""

Code: Select all

``````"
Created on Thu May 28 16:27:29 2020

@author: praf_malla@hotmail.com, Nepal
Corona virus time, Lockdown.
Stay home, Keep community safe.
"""
# =============================================================================
#  ANALYSIS OF RC COLUMN
# =============================================================================
#  Units in N, mm, sec
print("==========================")

from openseespy.postprocessing.Get_Rendering import *
from openseespy.opensees import *
wipe()

import numpy as np
import matplotlib.pyplot as plt
from math import asin, sqrt

print("Starting RC column analysis")
PI = 2 * asin(1.0)

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

# Set parameters for overall model geometry
height = 2000.0         # height of column

# Create nodes
node(1, 0.0, 0.0)
node(2, 0.0, height)

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

# Define cross-section for nonlinear columns
#  some parameters
cover    = 40   ;    # concrete cover
barDia   = 20   ;    # longitudinal bar dia
As = PI/4*barDia**2  # area of longitudinal bar

# Properties of concrete
# Define parameters of concrete
fc  = -30;             # Concrete compressive strength (+Tension, -compression)
Ec  = 5700*sqrt(-fc);  # Modulus of elasticity of concrete
e0  = 0.002;           # strain at compressive strength of concrete

# Define parameters of stirrups for confinement
fyh        = 415;    # yield strength of stirrups
diaStirrup = 8 ;     # Diameter of stirrups
(nLY, nLZ) = (2, 3)  # no of Legs of stirrup parallel to Y and Z - direction
sh         = 150 ;   # Spacing of stirrups

# Derived confinement parameters from stirrups
hP = colDepth-2*cover;  # Depth of concrete core to oustide of stirrups
bP = colWidth-2*cover;  # Width of concrete core to oustide of stirrups
Aso      = PI/4*diaStirrup**2;          # Area of stirrups
Vstirrup = Aso*(nLY*hP+nLZ*bP)*1000/sh; # Volume of stirrups per 1000 mm
Vcore    = hP*bP*1000;                  # Volume of concrete core in 1000 mm
Rhos     = Vstirrup/Vcore;              # ratio of volume of hoop steel to the volume of concrete core measured to the outside of hoops.

# Parameters of Confined concrete
Kfc   = 1+Rhos*fyh/(-fc);   # Ratio of confined to unconfined concrete strength (by Scott)
fc1C  = -Kfc*fc;            # Compressive strength of confined concrete
eps1C = Kfc*e0;             # Strain at confined compressive strength
fc2C  = 0.2*fc1C;           # stress of confined concrete at ultimate strain

e50uC = (0.0207+eps1C*fc1C)/(fc1C-6.89)
e50hC = 3/4*Rhos*sqrt(hP/sh)
Zc    = 0.5/(e50uC+e50hC-eps1C) ; # Descending branch of Confined concrete
eps2C = 0.8*fc1C/Zc+eps1C;        # strain at ultimate stress of confined concrete
ftUC  = 0.14*fc1C; EtsC = 2*ftUC/eps1C; lambdaC = 0.1

# Parameters of Unconfined concrete
fc1U  =-fc ;         # Compressive strength of unconfined concrete
eps1U = e0 ;         # strain at compressvie strength of unconfined concrete
fc2U  = 0.2*fc1U ;   # stress of unconfined concrete at ultimate strain

e50uU = (0.0207+eps1U*fc1U)/(fc1U-6.89)
e50hU = 0.0 ;                    # There is no sirrups i.e., Rhos = 0.0
Zu    = 0.5/(e50uU+e50hU-eps1U); # Descending branch of Unconfined concrete
eps2U = 0.8*fc1U/Zu+eps1U ;      # Strain at ultimate stress of unconfined concrete
ftUU  = 0.14*fc1U; EtsU = 2*ftUU/eps1U; lambdaU = 0.1

# ------------------------------------------
# CONCRETE                   tag
# Core concrete (confined)
uniaxialMaterial('Concrete02', 1, -fc1C, -eps1C, -fc2C, -eps2C, lambdaC , ftUC, EtsC)
# Cover concrete (unconfined)
uniaxialMaterial('Concrete02', 2, -fc1U, -eps1U, -fc2U, -eps2U, lambdaU , ftUU, EtsU)

print ("Parameters of Confined concrete", -fc1C, -eps1C, -fc2C, -eps2C)
print ("Parameters of UnConfined concrete", -fc1U, -eps1U, -fc2U, -eps2U)

# Plot stress-strain curve of Confined and Unconfined concrete
plt.figure()
sigC = (0, fc1C, fc2C); epsC = (0, eps1C, eps2C)
sigU = (0, fc1U, fc2U); epsU = (0, eps1U, eps2U)
plt.plot(epsC, sigC, linewidth=2.0)
plt.plot(epsU, sigU, linewidth=2.0)
plt.xlabel('strain')
plt.ylabel('compressive stress MPa')

# Longitudinal Reinforcing steel
fy = 400.0;     # Yield stress
E  = 200000.0;  # Young's modulus
Bs = 0.01;      # strain hardening ratio
(R0, cR1, cR2) = (18, 0.925, 0.15) # Steel parameters

# Variables derived for formation of section
coverY = colDepth / 2.0
coverZ = colWidth / 2.0
coreY  = coverY-cover
coreZ  = coverZ-cover
(nfCoreY, nfCoreZ, nfCoverY, nfCoverZ)= (20, 20, 20, 20)
secTag = 1
#                         tag
uniaxialMaterial('Steel02', 3, fy, E, Bs, R0, cR1, cR2)
section('Fiber', secTag)

# Create the concrete core fibers
patch('quad', 1, nfCoreZ, nfCoreY, -coreY, coreZ, -coreY, -coreZ, coreY, -coreZ, coreY, coreZ,)

# Create the concrete cover patches (right,left,bottom, top)
patch('quad', 2, 2, nfCoverY, -coverY, coverZ, - coreY, coreZ, coreY, coreZ, coverY, coverZ)
patch('quad', 2, 2, nfCoverY, -coreY, -coreZ, -coverY, -coverZ, coverY, -coverZ, coreY, -coreZ)
patch('quad', 2, nfCoverZ, 2, -coverY, coverZ, -coverY, -coverZ, -coreY, -coreZ, -coreY, coreZ)
patch('quad', 2, nfCoverZ, 2, coreY, coreZ, coreY, -coreZ, coverY, -coverZ, coverY, coverZ)

# Create the reinforcing fibers (top, middle, bottom)
layer('straight', 3, 3, As, coreY, coreZ, coreY, -coreZ)
layer('straight', 3, 2, As, 0.0, -coreZ, 0.0, coreZ)
layer('straight', 3, 3, As, -coreY, coreZ, -coreY, -coreZ)

# Define column elements
# ----------------------
geomTag = 1
# Geometry of column elements
#                tag
geomTransf('Corotational', geomTag)

np = 4 # Number of integration points along length of element
# Lobatto integration
#beamIntegration('Lobatto', 1, 1, np)

# Create the coulumns using Beam-column elements
#               e            tag ndI ndJ transfTag integrationTag
#eleType = 'forceBeamColumn'
#element(eleType, 1, 1, 2, 1, 1)

#       mixedBeamColumn tag iNode jNode numIntgrPts secTag transftag
eleType = 'mixedBeamColumn'
element(eleType, 1, 1, 2, np, secTag, geomTag, '-integration', 'Lobatto')

#  a parameter for the axial load
AxialLoad = 0.1*fc1U*colDepth*colWidth;  # 10% of axial capacity of columns

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

# Create nodal loads at nodes 3 & 4
#    nd  FX,  FY, MZ
# ------------------------------
# End of model generation
# ------------------------------

# ------------------------------
# Start of analysis generation
# ------------------------------
# Gravity Analysis
# Create the system of equation, a sparse solver with partial pivoting
Tol = 1.0e-6
system      ('BandGeneral')
constraints ('Plain')
numberer    ('RCM')
test        ('NormDispIncr', Tol, 10, 3)
algorithm   ('Newton')
analysis    ('Static')
analyze     (10)

# Display Model
plot_model()

print("==========================")
print("Gravity Analysis Completed")

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

print("Pushover Analysis Started!")
# Set the gravity loads to be constant & reset the time in the domain
# Set some parameters
# Set lateral load pattern with a Linear TimeSeries
pattern('Plain', 2, 1)

# Create nodal loads at nodes 3 & 4
#    nd    FX  FY  MZ

# ----------------------------------------------------
# Start of push over analysis
# ----------------------------------------------------

# Set some parameters
dU    = 1       # Displacement increment
dUMin = 1
dUMax = 1

(Tol, MaxIter) = (1.0e-12, 10)      # Tolerance and Maximum iteration
(controlNode, controlDOF) =(2, 1)   # Control node and control DOF
system       ('BandGeneral')
constraints  ('Plain')
numberer     ('RCM')
test         ('NormDispIncr', Tol, MaxIter)
analysis     ('Static');

# Change the integration scheme to be displacement control
#                             node dof init Jd min max
integrator   ('DisplacementControl', controlNode, controlDOF, dU, 1, dUMin,dUMax)
algorithm    ('Newton')

# Set some parameters
maxU        = 140.0  # Max displacement
currentDisp = nodeDisp(controlNode, controlDOF)
ok          = 0

# perform the analysis
Nsteps=int(maxU/dU)
import numpy as np
Nsteps=int(maxU/dU)             # Number of data to be stored
data=np.zeros( (Nsteps+1, 2) )  # Initial variable to store data
j = 0                           # Counter to store data
while (ok == 0 and currentDisp < maxU):
j+=1
analyze(1)
# if the analysis fails try initial tangent iteration
if ok != 0:
print("regular newton failed .. lets try an initail stiffness for this step")
test('NormDispIncr', 1.0e-6, 1000)
algorithm('ModifiedNewton', '-initial')
ok = analyze(1)
algorithm('Newton')

currentDisp = nodeDisp(controlNode, controlDOF)