Shear limit curve - Githubissues

Try it with the latest version of OpenSeesPy, then report back if there is an error or not.
Hi @mhscott @zhuminjie @u-anurag ,
I tried to use it, but it is still showing the warning "to be implemented". I am posting the Python code to reproduce the issue: (limit limitCurve('Shear') is available but limitCurve('Axial',...) is not available.
Can you point me to some resource on how to integrate an existing tcl function of OpneSees in OpenSeesPy
Note: This is the python version of the already existing Tcl code for the limit curve example. Please let me know if you need further input.
from openseespy.opensees import *
import math

# Number of analysis steps and step size
nSteps = 6000
dlamda = 0.1

# Loading option
loading = "cyclic"

# Analysis control option
control = "displ"

###########################
# Build model
###########################

wipe()
model('BasicBuilder', '-ndm', 2, '-ndf', 3)
print("HI")

################################
# Define nodal mesh and B.C.s
################################
L = 58.0

#    tag  X   Y
node(1, 0.0, 0.0)
node(2, 0.0, 0.0)
node(3, 0.0, L)
node(4, 0.0, L)
node(5, 0.0, L)

#   tag DX DY RZ
fix(1, 1, 1, 1)
fix(4, 0, 0, 1)
fix(5, 1, 1, 1)

##############################
# Create column section
##############################

# source('CenterColSecFiber.tcl')
# Parameters
h = 9.0
b = h
beta1 = 100.0
beta2 = 100.0
spall = "nospalling"
alphaS = 0.015
nfCore = 28
nfCover = 4
c = 1.0
z = h / 2.0 - c
y = b / 2.0 - c
nz = -z
ny = -y
zc = z + c
yc = y + c
nzc = -zc
nyc = -yc
fc = -3.517
fc = -3.517 * beta2 / 100.0
fcu = fc * beta1 / 100.0
fy = 69.5
Es = 29000.0
esy = fy / Es
esu = 0.15
fu = fy + alphaS * Es * (esu - esy)
Ec = 3400
A = h * b
L = 58
Ig = h * h * h * b / 12
Mcr = 228.0
Kcr = 1.2e-4
My = 554.0
Ky = 6.3e-4
Ry = Ky * L / 6.0
EI = My / Ky
alpha = 0.05
Ku = 0.05
Mu = My + alpha * EI * (Ku - Ky)
Ru = 0.5
pinchX = 0.5
pinchY = 0.4
damage1 = 0.0
damage2 = 0.0
beta = 0.4
Ic = EI / Ec
slipStiff = 91000
Rslipy = My / slipStiff
alphaSlip = alpha / (1.0 + (1.0 - alpha) * (slipStiff * L / 6 / EI))
Rslipu = Rslipy + (Mu - My) / (alphaSlip * slipStiff)

# Create axial failure spring
# source('CenterColAxialSpring.tcl')
#axial spring

# Parameters
Fsw = 11.87 / 1.0
axialElasticSlope = 99.0 * Ec * A / L  # 99 times more rigid than column
axialNegSlope = -90.0
Pr = 5.0
P1 = 65.0
P2 = 75.0
P3 = 85.0

# Define limit curve for axial spring
axialCurveTag = 2
bcTag = 99

# Define LimitStateMaterial
axialFailTag = 8

# Define the limit curve
limit_curve_params = [Fsw, axialNegSlope, Pr, 2, 2, 1, 4, 1, 2, 0.0, 0]
limitCurve('Axial', axialCurveTag, bcTag, *limit_curve_params)

# Convert axial loads used for setting initial elastic slope
P1 /= axialElasticSlope
P2 /= axialElasticSlope
P3 /= axialElasticSlope

# Define LimitStateMaterial
uniaxialMaterial('LimitState', axialFailTag, P1, P1, P2, P2, P3, P3, -P1, -P1, -P2, -P2, -P3, -P3, 0.5, 0.5, 0.0, 0.0, 0.0, axialCurveTag, 1)

# Create shear failure spring 
# source('CenterColShearSpring.tcl')
# Parameters
#shear spring
rigidSlope = 1700
negSlope = -8
Vr = 3.0
Vi1 = 25.0
Vi2 = 30.0
Vi3 = 45.0
kf = 24.7

# Define limit curve using shear drift model
shearCurveTag = 1
bcTag = 99

# Define HystereticMaterial
shearTag = 4

# Define the limit curve
limit_curve_params = [0.0018, 3517.0, 9.0, 9.0, 7.75, 11.87, kf, Vr, 2, 0, 1, 4, 1, 2, 0.0]
limitCurve('Shear', shearCurveTag, bcTag, *limit_curve_params)

# Convert strengths for initial response
Vi1 /= rigidSlope
Vi2 /= rigidSlope
Vi3 /= rigidSlope

# Define HystereticMaterial
uniaxialMaterial('LimitState', shearTag, Vi1, Vi1, Vi2, Vi2, Vi3, Vi3, -Vi1, -Vi1, -Vi2, -Vi2, -Vi3, -Vi3, 0, 0, 0, 0, 0, shearCurveTag, 2, 0)

section('Aggregator', shearAxialOnlySec, shearTag, 'Vy', axialFailTag, 'P')

# Create fiber section
# Define uniaxialMaterials
#                           tag         f'c    epsc   f'cu   epscu
uniaxialMaterial('Concrete01', coreTag, fc, -0.002, fcu, -0.0052)

if spall == "spalling":
    uniaxialMaterial('Concrete01', coverTag, fc, -0.002, 0.0, -0.0060)
elif spall == "nospalling":
    uniaxialMaterial('Concrete01', coverTag, fc, -0.002, fcu, -0.0052)
else:
    print("Invalid spalling option: {}".format(spall))

#                           tag          fy     E      hardening ratio
# uniaxialMaterial Steel02     $steelTag   69.5   29000  $alphaS
uniaxialMaterial('Hysteretic', steelTag, fy, esy, fu, esu, -fy, -esy, -fu, -esu, 1.0, 1.0, 0, 0)

# Define the fiber section
section('Fiber', flexSec)
# Define the concrete patch with fibers for unidirectional bending
patch('quadr', coreTag, 1, nfCore, ny, z, ny, nz, y, nz, y, z)
# Define the four cover patches
patch('quadr', coverTag, 1, nfCover, nyc, zc, nyc, nzc, ny, nz, ny, z)
patch('quadr', coverTag, 1, nfCover, y, z, y, nz, yc, nzc, yc, zc)
patch('quadr', coverTag, 1, nfCore, ny, nz, nyc, nzc, yc, nzc, y, nz)
patch('quadr', coverTag, 1, nfCore, nyc, zc, ny, z, y, z, yc, zc)
# Define the reinforcement explicitly using fiber command
#     yloc  zloc  area matTag
fiber(-3.250, 3.250, 0.2, steelTag, 0)
fiber(-3.250, -3.250, 0.2, steelTag, 0)
fiber(3.250, -3.250, 0.2, steelTag, 0)
fiber(3.250, 3.250, 0.2, steelTag, 0)
fiber(-3.187, 0.0, 0.31, steelTag, 0)
fiber(3.187, 0.0, 0.31, steelTag, 0)
fiber(0.0, -3.187, 0.31, steelTag, 0)
fiber(0.0, 3.187, 0.31, steelTag, 0)

Acenter = A
Eccenter = Ec
Iccenter = Ic

#################################
# Define the beam-column element
#################################
geomTransf('PDelta', 1)
nint = 5
element('nonlinearBeamColumn', bcTag, 2, 3, nint, flexSec, 1, '-iter', 5, '1e-15')

####################################
# Define the zero-length end springs
####################################
# rigid material
uniaxialMaterial('Elastic', rigidMatTag, 9.9e9)

# bottom of column slip spring 
element('zeroLength', 1, 1, 2, '-mat', rigidMatTag, rigidMatTag, centerSlipTag, '-dir', 1, 2, 6)

# top of column springs incl shear and axial
element('zeroLength', 3, 3, 4, '-mat', shearTag, axialFailTag, centerSlipTag, '-dir', 1, 2, 6)

# soft spring to take axial load after failure
uniaxialMaterial('Elastic', softMatTag, 99.9)
element('zeroLength', 4, 4, 5, '-mat', softMatTag, '-dir', 2)

# Recorders
# record displacements for node with applied displacement
recorder('Node', 'node4basic.out', 'disp', '-time', '-node', 4, '-dof', 1, 2, 3)

# record end section forces and deformations
recorder('Element', bcTag, '-time', '-file', 'secforce1basic.out', 'section', 1, 'force')
recorder('Element', bcTag, '-time', '-file', 'secdeform1basic.out', 'section', 1, 'deformation')

# record beam-column element forces and rotations
recorder('Element', 4, '-time', '-file', 'axialSpr.out', 'force')
recorder('Element', bcTag, '-time', '-file', 'eleforcebasic.out', 'force')
recorder('Element', bcTag, '-time', '-file', 'elerotbasic.out', 'rotation')

# Apply gravity loads
# Initial axial load 
P = -70.0
# Constant gravity load
pattern('Plain', 1, 'Constant', {
    load(4, 0.0, P, 0.0)
})
# Analysis options for gravity loads
initialize()
system('ProfileSPD')
constraints('Penalty', 1.0e12, 1.0e12)
numberer('RCM')
test('NormDispIncr', 1.0e-8, 25, 0)
algorithm('Newton')
integrator('LoadControl', 0, 1, 0, 0)
analysis('Static')
# Apply gravity load in one step
analyze(1)
wipeAnalysis()
# Apply transverse displacements
if loading == "push":
    # Define displacements of node 4 for pushover
    du = 0.01
    pattern('Plain', 2, f"Linear -factor {du}", {
        sp(4, 1, 1.0)
    })
elif loading == "cyclic":
    path = [1.2557524028480316e-003, 1.2557524028480316e-003, 3.0608252487525078e-003, 1.5641001656305775e-003,
            1.1015791633965932e-003, -1.8313432524451514e+000]
    pattern('Plain', 2, f"Series -dt 1.0 -values {path}", {
        sp(4, 1, 1.0)
    })
else:
    print(f"Invalid loading option: {loading}")
# Set Analysis options for transverse loading
# Create the system of equation
system('BandGeneral')
# Create the DOF numberer, the reverse Cuthill-McKee algorithm
numberer('RCM')
# Create the integration scheme
if control == "load":
    # LoadControl scheme using constant steps of dlamda
    integrator('LoadControl', dlamda, 1, dlamda, dlamda)
    constraints('Penalty', 1.0e14, 1.0e14)
elif control == "displ":
    # DisplacementControl scheme with constant loading rate
    integrator('DisplacementControl', 2, 1, du, 1, du, du)
    constraints('Penalty', 1.0e14, 1.0e14)
else:
    print(f"Invalid control option: {control}")
# Create the convergence test, Norm of the displacement increment
test('NormDispIncr', 1.0e-10, 25, 0)
# Create the solution algorithm, a Newton-Rahpson algorithm is created
algorithm('Newton')
# Create the analysis object
analysis('Static')
# Analyze model with transverse loading
ok = 0
n = 1
while n < nSteps and ok == 0:
    ok = analyze(1)
    if ok != 0:
        test('NormDispIncr', 1.0e-10, 100000, 1)  # increase max number of iterations
        algorithm('ModifiedNewton', '-initial')  # use initial elastic stiffness for NewtonRaphson
        # print(f"Time Newton Initial {getTime()}")  # output time algorithm was changed
        ok = analyze(1)  # analyze 1 step with new settings
        # ans = input()  # pauses tcl script
        algorithm('Newton')  # restore algorithm to Newton with current stiffness
        test('NormDispIncr', 1.0e-10, 25, 0)  # restore max number of iterations
    n += 1

if ok != 0:
    print("ANALYSIS FAILED")
else:
    print("SUCCESSFUL")
wipe()
zhuminjie / OpenSeesPyDoc

Shear limit curve #238
