The OpenSees Community

Yes I had tried each line separately and in combination to check the node reactions. There is a large increase in reaction on application of the load, which then dissipates, I guess as the structure reacts and displaces. But in a simpler model, the affect would be much clearer. Many thanks again for your response!

Great - thankyou!!!

Hi All,
In parallel computing I am receiving the following error:

WARNING MumpsParallelSolver::solve(void)- Error -9 returned in substitution dmumps()
Work array too small; use -ICNTL14 option, the default is -ICNTL 20 make 20 larger

However if I try to increase the ICNTL14 in openseepy by using e.g.:
ops.system('Mumps','-ICNTL14', 80)

But still get the same error message, unless I reduce the number of processors.

Does anyone have an example of calling this in OpenseePy?

Many thanks in advance.
Kind regards

Hi All,
I am receiving the following error:

WARNING MumpsParallelSolver::solve(void)- Error -9 returned in substitution dmumps()
Work array too small; use -ICNTL14 option, the default is -ICNTL 20 make 20 larger

However if I try to increase the ICNTL14 in openseepy by using e.g.:
ops.system('Mumps','-ICNTL14', 80)

But still get the same error message, unless I reduce the number of processors.

Does anyone have an example of calling this in OpenseePy?

Many thanks in advance.
Kind regards

Hi All,
Does anyone have experience of applying simultaneous multidirectional loads in a transient analysis to a structure in Opensees?
I would like to apply loading at 3 points onto the structure, and have time series for the x,y, and z component of the loads.

To do this I am using the below in my script, these seem to be ok if you print the model (see print out below the script), but I want to check if all these load patterns are actually being applied simultaneously, or it is just the last load pattern before the analysis call that is actually applied.
Many thanks if you have tried this or similar / have any suggestions.

# Time series
ops.timeSeries('Path', 1000, '-values','-time', '-filePath', 'T1_x_Int_001.txt', '-fileTime', 'time_0.001.txt')
ops.timeSeries('Path', 2000, '-values','-time', '-filePath', 'T1_y_Int_001.txt', '-fileTime', 'time_0.001.txt')
ops.timeSeries('Path', 3000, '-values','-time', '-filePath', 'T1_z_Int_001.txt', '-fileTime', 'time_0.001.txt')

ops.timeSeries('Path', 4000, '-values','-time', '-filePath', 'T2_x_Int_001.txt', '-fileTime', 'time_0.001.txt')
ops.timeSeries('Path', 5000, '-values','-time', '-filePath', 'T2_y_Int_001.txt', '-fileTime', 'time_0.001.txt')
ops.timeSeries('Path', 6000, '-values','-time', '-filePath', 'T2_z_Int_001.txt', '-fileTime', 'time_0.001.txt')

ops.timeSeries('Path', 7000, '-values','-time', '-filePath', 'T3_x_Int_001.txt', '-fileTime', 'time_0.001.txt')
ops.timeSeries('Path', 8000, '-values','-time', '-filePath', 'T3_y_Int_001.txt', '-fileTime', 'time_0.001.txt')
ops.timeSeries('Path', 9000, '-values','-time', '-filePath', 'T3_z_Int_001.txt', '-fileTime', 'time_0.001.txt')

#T1

ops.pattern('Plain', 1, 1000)
ops.load(31755, 1.0, 0.0, 0.0) # Load in X direction

ops.pattern('Plain', 2, 2000)
ops.load(31755, 0.0, 1.0, 0.0) # Load in Y direction

ops.pattern('Plain', 3, 3000)
ops.load(31755, 0.0, 0.0, 1.0) # Load in Z direction


#T2

ops.pattern('Plain', 4, 4000)
ops.load(32045, 1.0, 0.0, 0.0) # Load in X direction

ops.pattern('Plain', 5, 5000)
ops.load(32045, 0.0, 1.0, 0.0) # Load in Y direction

ops.pattern('Plain', 6, 6000)
ops.load(32045, 0.0, 0.0, 1.0) # Load in Z direction


#T3

ops.pattern('Plain', 7, 7000)
ops.load(31947, 1.0, 0.0, 0.0) # Load in X direction

ops.pattern('Plain', 8, 8000)
ops.load(31947, 0.0, 1.0, 0.0) # Load in Y direction

ops.pattern('Plain', 9, 9000)
ops.load(31947, 0.0, 0.0, 1.0) # Load in Z direction

ops.analyze(10000, 0.001)

Model print out:

LOAD PATTERNS: numPatterns: 9

numComponents: 9
Load Pattern: 1
Scale Factor: 1
Path Time Series: constant factor: 1 Nodal Loads:
numComponents: 1
Nodal Load: 31755 load : 1 0 0
Elemental Loads:
numComponents: 0
Single Point Constraints:
numComponents: 0

Load Pattern: 2
Scale Factor: 1
Path Time Series: constant factor: 1 Nodal Loads:
numComponents: 1
Nodal Load: 31755 load : 0 1 0
Elemental Loads:
numComponents: 0
Single Point Constraints:
numComponents: 0

Load Pattern: 3
Scale Factor: 1
Path Time Series: constant factor: 1 Nodal Loads:
numComponents: 1
Nodal Load: 31755 load : 0 0 1
Elemental Loads:
numComponents: 0
Single Point Constraints:
numComponents: 0

Load Pattern: 4
Scale Factor: 1
Path Time Series: constant factor: 1 Nodal Loads:
numComponents: 1
Nodal Load: 32045 load : 1 0 0
Elemental Loads:
numComponents: 0
Single Point Constraints:
numComponents: 0

Load Pattern: 5
Scale Factor: 1
Path Time Series: constant factor: 1 Nodal Loads:
numComponents: 1
Nodal Load: 32045 load : 0 1 0
Elemental Loads:
numComponents: 0
Single Point Constraints:
numComponents: 0

Load Pattern: 6
Scale Factor: 1
Path Time Series: constant factor: 1 Nodal Loads:
numComponents: 1
Nodal Load: 32045 load : 0 0 1
Elemental Loads:
numComponents: 0
Single Point Constraints:
numComponents: 0

Load Pattern: 7
Scale Factor: 1
Path Time Series: constant factor: 1 Nodal Loads:
numComponents: 1
Nodal Load: 31947 load : 1 0 0
Elemental Loads:
numComponents: 0
Single Point Constraints:
numComponents: 0

Load Pattern: 8
Scale Factor: 1
Path Time Series: constant factor: 1 Nodal Loads:
numComponents: 1
Nodal Load: 31947 load : 0 1 0
Elemental Loads:
numComponents: 0
Single Point Constraints:
numComponents: 0

Load Pattern: 9
Scale Factor: 1
Path Time Series: constant factor: 1 Nodal Loads:
numComponents: 1
Nodal Load: 31947 load : 0 0 1
Elemental Loads:
numComponents: 0
Single Point Constraints:
numComponents: 0
PARAMETERS: numParameters: 0
numComponents: 0

If anyone is wondering about this you can add to the model elements with different dof, by calling the model command, see text within:

https://opensees.github.io/OpenSeesDocu ... model.html

Also, you can then join elements together with different dof by using the Equal dof:

https://openseespydoc.readthedocs.io/en ... alDOF.html

A geotech example is shown here:

https://opensees.berkeley.edu/wiki/inde ... ifuge_Test

Hi All,

Does anyone have any suggestions for how to use SSPBrick (3dof) and SSPBrickUP (3dof +pwp dof) elements in the same model?

I would like to use SSPBrick to model the pile and SSPBrickUP for the soil.

I have tried a few different dof setting in the model set up etc but none have been successful.

Any suggestions welcome,

Many thanks!

Hi All,

Has anyone (recently) used SAniSandMS successfully in opensees/ openseepy?

The below script is for a single element drained triaxial. It works OK for the older version of Sanisand04 (Manzari Dafalias model), but will not work with SAniSandMS. There is sometimes an error:

SAniSandMS::setResponse -- Unrecognized response option "fabric"

In addition to the above (which could certainly be my modelling error!), if you try the tcl example for SAniSandMS provided in opensees (3.5) – copy and paste it, the output does not resemble the suggested visual output. The cyclic loading is applied properly, but the principle strain generation does not accumulate (ratchet), but cycles at a strain increment of around -4e-3. This is different from the ratcheting displayed in the output of the example and expected performance of the model. This tcl example is located:

https://opensees.github.io/OpenSeesDocu ... andMS.html

So I would be interested to hear if anyone has successfully used SAniSandMS in Openseepy recently?

And, if they have any suggestions ( I am trying very small time steps) …. or even perhaps there is something wrong with the model at the moment/ function parameter list (fabric deprecated?), as in the tcl example, as far I can see the output does not seem to produce the results shown in the example.
Many thanks in advance!

The drained TXL script for Manzari Dafalias (working) / SAniSandMS (strange results) is below (NB long running time due to small dT):


import openseespy.opensees as op
import time
import numpy as np
import matplotlib.pyplot as plt


# Test Specific parameters
pConf = -300.0
devDisp = -0.3
perm = 1.0e-10
vR = 0.8

damp = 0.1
omega1 = 0.0157
omega2 = 64.123
a1 = 2.0 * damp / (omega1 + omega2)
a0 = a1 * omega1 * omega2

op.wipe()

# Create a 3D model with 4 Degrees of Freedom

op.model('BasicBuilder', '-ndm', 3, '-ndf', 4)

# Create nodes

op.node(1, 1.0, 0.0, 0.0)
op.node(2, 1.0, 1.0, 0.0)
op.node(3, 0.0, 1.0, 0.0)
op.node(4, 0.0, 0.0, 0.0)
op.node(5, 1.0, 0.0, 1.0)
op.node(6, 1.0, 1.0, 1.0)
op.node(7, 0.0, 1.0, 1.0)
op.node(8, 0.0, 0.0, 1.0)

# Create Fixities (dof 4 = 1 = drained)

op.fix (1, 0, 1, 1, 1)
op.fix (2, 0, 0, 1, 1)
op.fix (3, 1, 0, 1, 1)
op.fix (4, 1, 1, 1, 1)
op.fix (5, 0, 1, 0, 1)
op.fix (6, 0, 0, 0, 1)
op.fix (7, 1, 0, 0, 1)
op.fix (8, 1, 1, 0, 1)

# Create material

G0 = 110.0 # [Adimensional]
nu = 0.05 # [Adimensional]
e_init = 0.72 # [Adimensional]
Mc = 1.27 # [Adimensional]
c = 0.712 # [Adimensional]
lambda_c = 0.049 # [Adimensional]
e0 = 0.845 # [Adimensional]
ksi = 0.27 # [Adimensional]
P_atm = 101.3 # [kPa]
m = 0.01 # [Adimensional]
h0 = 5.95 # [Adimensional]
ch = 1.01 # [Adimensional]
nb = 2.0 # [Adimensional]
A0 = 1.06 # [Adimensional]
nd = 1.17 # [Adimensional]
z_max = 4 # For SAniSand [Adimensional]
cz = 0 # For SAniSand [Adimensional]
mu0 = 260.0 # For SAniSand [Adimensional]
zeta = 0.0005 # For SAniSand [Adimensional]
beta = 1.0 # For SAniSand [Adimensional]
w1 = 0.5
w2 = 2
#fabric_flag
#flow_flag
Den = 1.584 # [Mg/m^3]
intScheme = 3 # Corresponds to Modified-Euler integration scheme
TanType = 1 # 0: elastic stiffness, 1: continuum elastoplastic stiffness
JacoType = 1 # Not used in explicit methods
TolF = 1.0e-6 # Tolerances, not used in explicit
TolR = 1.0e-6 # Tolerances, not used in explicit



#Reference atmospheric pressure
P_ref = -101.3

# Create material
#nDMaterial SAniSandMS $matTag $G0 $nu $e_init $Mc $c $lambda_c $e0 $ksi $P_atm $m $h0 $ch $nb $A0 $nd $zeta $mu0 $beta $Den $fabric_flag $flow_flag $intScheme $TanType $JacoType $TolF $TolR

#nDMaterial SAniSandMS 1 $G0 $nu $e_init $Mc $c $lambda_c $e0 $ksi $P_atm $m $h0 $ch $nb $A0 $nd $zeta $mu0 $beta $Den $intScheme $TanType $JacoType $TolF $TolR

op.nDMaterial('SAniSandMS', 1, G0, nu, e_init, Mc, c, lambda_c, e0, ksi, P_atm, m, h0, ch, nb, A0, nd, zeta, mu0, beta, Den, intScheme, TanType, JacoType, TolF, TolR)
op.type = "RK"

# nDmaterial ManzariDafalias $matTag $G0 $nu $e_init $Mc $c $lambda_c $e0 $ksi $P_atm $m $h0 $ch $nb $A0 $nd $z_max $cz $Den
#op.nDMaterial('ManzariDafalias', 1, 125, 0.05, vR, 1.25, 0.712, 0.019, 0.934, 0.7, 100, 0.01, 7.05, 0.968, 1.1, 0.704, 3.5, 4, 600, 1.42)


# Create element
op.element('SSPbrickUP', 1, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2.2e6, 1.0, perm, perm, perm, vR, 1.5e-9)

# Create recorders
op.recorder('Node', '-file', 'disp.txt', '-time', '-nodeRange', 1, 8, '-dof', 1, 2, 3, 'disp')
op.recorder('Node', '-file', 'press.txt', '-time', '-nodeRange', 1, 8, '-dof', 4, 'vel')

op.recorder('Element','-file', 'Strain.txt','-time','-ele', 1, 'strain')
op.recorder('Element','-file', 'Stress.txt','-time','-ele', 1, 'stress')
op.recorder('Element','-file', 'Alpha.txt','-time','-ele', 1, 'alpha')
op.recorder('Element','-file', 'Fabric.txt','-time','-ele', 1, 'fabric')


# Create analysis

op.constraints('Transformation')
op.test('NormDispIncr', 1.0e-5, 20, 0)
op.algorithm('Newton')
op.numberer('RCM')
op.system('BandGeneral')
op.integrator('Newmark', 0.5, 0.25)
op.analysis('Transient')

# Apply confinement pressure
pNode = pConf / 4.0

op.timeSeries('Path', 1, '-values', 0, 1, 1, '-time', 0.0, 1000.0, 10000000000.)
op.pattern("Plain", 1, 1, '-factor',1.0)

op.load (1, pNode, 0.0, 0.0, 0.0),
op.load (2, pNode, pNode, 0.0, 0.0),
op.load (3, 0.0, pNode, 0.0, 0.0),
op.load (4, 0.0, 0.0, 0.0, 0.0),
op.load (5, pNode, 0.0, pNode, 0.0),
op.load (6, pNode, pNode, pNode, 0.0),
op.load (7, 0.0, pNode, pNode, 0.0),
op.load (8, 0.0, 0.0, pNode, 0.0)

# attempt at psuedo static loading
op.analyze(10000,0.1)

# Read vertical displacement of top plane
vertDisp = op.nodeDisp(5, 3)

# Apply deviatoric strain

lValues = [1, 1 + devDisp/vertDisp, 1 + devDisp/vertDisp]

# Define the time series
op.timeSeries('Path', 2, '-values', *lValues, '-time', 1000.0, 2000.0, 10000000, '-factor', 1.0)

# Create the loading pattern
op.pattern('Plain', 2, 2, '-factor', 1)

op.sp(5, 3, vertDisp)
op.sp(6, 3, vertDisp)
op.sp(7, 3, vertDisp)
op.sp(8, 3, vertDisp)

# Set number and length of (pseudo)time steps

dT = 0.001
numStep = 1000000

# Initialize the remaining steps and success flag
remStep = numStep
success = 0

def subStepAnalyze(dT, subStep):
if subStep > 10:
return -10

for i in range(1, 3):
print(f"Try dT = {dT}")
success = op.analyze(1, dT)
if success != 0:
success = subStepAnalyze(dT/2.0, subStep + 1)
if success == -10:
print("Did not converge.")
return success
else:
if i == 1:
print(f"Substep {subStep} : Left side converged with dT = {dT}")
else:
print(f"Substep {subStep} : Right side converged with dT = {dT}")

return success

import time

print("Start analysis")
startT = time.time()

while success != -10:
subStep = 0
success = op.analyze(remStep, dT)
if success == 0:
print("Analysis Finished")
break
else:
curTime = op.getTime()
print(f"Analysis failed at {curTime}. Try substepping.")
success = subStepAnalyze(dT / 2.0, subStep + 1)
curStep = int((curTime - 20000) / dT + 1)
remStep = int(numStep - curStep)
print(f"Current step: {curStep}, Remaining steps: {remStep}")

endT = time.time()
print(f"Loading analysis execution time: {endT - startT} seconds.")

op.wipe()

#%%PostProcessing
import pandas as pd
df_stress = pd.read_csv('stress.txt', sep=" ", header=None)
df_strain = pd.read_csv('strain.txt', sep=" ", header=None)
#df_PP = pd.read_csv('CycPP.txt', sep=" ", header=None)

#stress yy = stress xx
Stress_yy = df_stress.iloc[:, 2].to_numpy()*(-1.0) #compression is positive
Principle_Stress = df_stress.iloc[:, 3].to_numpy()*-1.0
Principle_Strain = df_strain.iloc[:, 3].to_numpy()*-1.0

fig = plt.figure(figsize=(10,18))
ax0 = fig.add_subplot(211)
ax0.plot(Principle_Strain, Principle_Stress, label='stress-strain', linewidth=1.2, color='red')
ax0.set_xlabel("Principle strain,%", fontsize=24)
ax0.set_ylabel("Principle stress, kPa", fontsize=24)
ax0.legend(fontsize=24)

ax1 = fig.add_subplot(212)
ax1.plot(Stress_yy, Principle_Stress, label='Stress Path', linewidth=1.2, color='red')
ax1.set_xlabel("Stress yy = xx, kPa", fontsize=24)
ax1.set_ylabel("Principle stress, kPa", fontsize=24)
ax1.legend(fontsize=24)
plt.show()
plt.close()

Hi All,

Has anyone (recently) used SAniSandMS successfully in opensees/ openseepy?

The below script is for a single element drained triaxial. It works OK for the older version of Sanisand04 (Manzari Dafalias model), but will not work with SAniSandMS. There is sometimes an error:

SAniSandMS::setResponse -- Unrecognized response option "fabric"

In addition to the above (which could certainly be my modelling error!), if you try the tcl example for SAniSandMS provided in opensees (3.5) – copy and paste it, the output does not resemble the suggested visual output. The cyclic loading is applied properly, but the principle strain generation does not accumulate (ratchet), but cycles at a strain increment of around -4e-3. This is different from the ratcheting displayed in the output of the example and expected performance of the model. This tcl example is located:

https://opensees.github.io/OpenSeesDocu ... andMS.html

So I would be interested to hear if anyone has successfully used SAniSandMS in Openseepy recently?

And, if they have any suggestions ( I am trying very small time steps) …. or even perhaps there is something wrong with the model at the moment/ function parameter list (fabric deprecated?), as in the tcl example, as far I can see the output does not seem to produce the results shown in the example.
Many thanks in advance!

The drained TXL script for Manzari Dafalias (working) / SAniSandMS (strange results) is below (NB long running time due to small dT):


import openseespy.opensees as op
import time
import numpy as np
import matplotlib.pyplot as plt


# Test Specific parameters
pConf = -300.0
devDisp = -0.3
perm = 1.0e-10
vR = 0.8

damp = 0.1
omega1 = 0.0157
omega2 = 64.123
a1 = 2.0 * damp / (omega1 + omega2)
a0 = a1 * omega1 * omega2

op.wipe()

# Create a 3D model with 4 Degrees of Freedom

op.model('BasicBuilder', '-ndm', 3, '-ndf', 4)

# Create nodes

op.node(1, 1.0, 0.0, 0.0)
op.node(2, 1.0, 1.0, 0.0)
op.node(3, 0.0, 1.0, 0.0)
op.node(4, 0.0, 0.0, 0.0)
op.node(5, 1.0, 0.0, 1.0)
op.node(6, 1.0, 1.0, 1.0)
op.node(7, 0.0, 1.0, 1.0)
op.node(8, 0.0, 0.0, 1.0)

# Create Fixities (dof 4 = 1 = drained)

op.fix (1, 0, 1, 1, 1)
op.fix (2, 0, 0, 1, 1)
op.fix (3, 1, 0, 1, 1)
op.fix (4, 1, 1, 1, 1)
op.fix (5, 0, 1, 0, 1)
op.fix (6, 0, 0, 0, 1)
op.fix (7, 1, 0, 0, 1)
op.fix (8, 1, 1, 0, 1)

# Create material

G0 = 110.0 # [Adimensional]
nu = 0.05 # [Adimensional]
e_init = 0.72 # [Adimensional]
Mc = 1.27 # [Adimensional]
c = 0.712 # [Adimensional]
lambda_c = 0.049 # [Adimensional]
e0 = 0.845 # [Adimensional]
ksi = 0.27 # [Adimensional]
P_atm = 101.3 # [kPa]
m = 0.01 # [Adimensional]
h0 = 5.95 # [Adimensional]
ch = 1.01 # [Adimensional]
nb = 2.0 # [Adimensional]
A0 = 1.06 # [Adimensional]
nd = 1.17 # [Adimensional]
z_max = 4 # For SAniSand [Adimensional]
cz = 0 # For SAniSand [Adimensional]
mu0 = 260.0 # For SAniSand [Adimensional]
zeta = 0.0005 # For SAniSand [Adimensional]
beta = 1.0 # For SAniSand [Adimensional]
w1 = 0.5
w2 = 2
#fabric_flag
#flow_flag
Den = 1.584 # [Mg/m^3]
intScheme = 3 # Corresponds to Modified-Euler integration scheme
TanType = 1 # 0: elastic stiffness, 1: continuum elastoplastic stiffness
JacoType = 1 # Not used in explicit methods
TolF = 1.0e-6 # Tolerances, not used in explicit
TolR = 1.0e-6 # Tolerances, not used in explicit



#Reference atmospheric pressure
P_ref = -101.3

# Create material
#nDMaterial SAniSandMS $matTag $G0 $nu $e_init $Mc $c $lambda_c $e0 $ksi $P_atm $m $h0 $ch $nb $A0 $nd $zeta $mu0 $beta $Den $fabric_flag $flow_flag $intScheme $TanType $JacoType $TolF $TolR

#nDMaterial SAniSandMS 1 $G0 $nu $e_init $Mc $c $lambda_c $e0 $ksi $P_atm $m $h0 $ch $nb $A0 $nd $zeta $mu0 $beta $Den $intScheme $TanType $JacoType $TolF $TolR

op.nDMaterial('SAniSandMS', 1, G0, nu, e_init, Mc, c, lambda_c, e0, ksi, P_atm, m, h0, ch, nb, A0, nd, zeta, mu0, beta, Den, intScheme, TanType, JacoType, TolF, TolR)
op.type = "RK"

# nDmaterial ManzariDafalias $matTag $G0 $nu $e_init $Mc $c $lambda_c $e0 $ksi $P_atm $m $h0 $ch $nb $A0 $nd $z_max $cz $Den
#op.nDMaterial('ManzariDafalias', 1, 125, 0.05, vR, 1.25, 0.712, 0.019, 0.934, 0.7, 100, 0.01, 7.05, 0.968, 1.1, 0.704, 3.5, 4, 600, 1.42)


# Create element
op.element('SSPbrickUP', 1, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2.2e6, 1.0, perm, perm, perm, vR, 1.5e-9)

# Create recorders
op.recorder('Node', '-file', 'disp.txt', '-time', '-nodeRange', 1, 8, '-dof', 1, 2, 3, 'disp')
op.recorder('Node', '-file', 'press.txt', '-time', '-nodeRange', 1, 8, '-dof', 4, 'vel')

op.recorder('Element','-file', 'Strain.txt','-time','-ele', 1, 'strain')
op.recorder('Element','-file', 'Stress.txt','-time','-ele', 1, 'stress')
op.recorder('Element','-file', 'Alpha.txt','-time','-ele', 1, 'alpha')
op.recorder('Element','-file', 'Fabric.txt','-time','-ele', 1, 'fabric')


# Create analysis

op.constraints('Transformation')
op.test('NormDispIncr', 1.0e-5, 20, 0)
op.algorithm('Newton')
op.numberer('RCM')
op.system('BandGeneral')
op.integrator('Newmark', 0.5, 0.25)
op.analysis('Transient')

# Apply confinement pressure
pNode = pConf / 4.0

op.timeSeries('Path', 1, '-values', 0, 1, 1, '-time', 0.0, 1000.0, 10000000000.)
op.pattern("Plain", 1, 1, '-factor',1.0)

op.load (1, pNode, 0.0, 0.0, 0.0),
op.load (2, pNode, pNode, 0.0, 0.0),
op.load (3, 0.0, pNode, 0.0, 0.0),
op.load (4, 0.0, 0.0, 0.0, 0.0),
op.load (5, pNode, 0.0, pNode, 0.0),
op.load (6, pNode, pNode, pNode, 0.0),
op.load (7, 0.0, pNode, pNode, 0.0),
op.load (8, 0.0, 0.0, pNode, 0.0)

# attempt at psuedo static loading
op.analyze(10000,0.1)

# Read vertical displacement of top plane
vertDisp = op.nodeDisp(5, 3)

# Apply deviatoric strain

lValues = [1, 1 + devDisp/vertDisp, 1 + devDisp/vertDisp]

# Define the time series
op.timeSeries('Path', 2, '-values', *lValues, '-time', 1000.0, 2000.0, 10000000, '-factor', 1.0)

# Create the loading pattern
op.pattern('Plain', 2, 2, '-factor', 1)

op.sp(5, 3, vertDisp)
op.sp(6, 3, vertDisp)
op.sp(7, 3, vertDisp)
op.sp(8, 3, vertDisp)

# Set number and length of (pseudo)time steps

dT = 0.001
numStep = 1000000

# Initialize the remaining steps and success flag
remStep = numStep
success = 0

def subStepAnalyze(dT, subStep):
if subStep > 10:
return -10

for i in range(1, 3):
print(f"Try dT = {dT}")
success = op.analyze(1, dT)
if success != 0:
success = subStepAnalyze(dT/2.0, subStep + 1)
if success == -10:
print("Did not converge.")
return success
else:
if i == 1:
print(f"Substep {subStep} : Left side converged with dT = {dT}")
else:
print(f"Substep {subStep} : Right side converged with dT = {dT}")

return success

import time

print("Start analysis")
startT = time.time()

while success != -10:
subStep = 0
success = op.analyze(remStep, dT)
if success == 0:
print("Analysis Finished")
break
else:
curTime = op.getTime()
print(f"Analysis failed at {curTime}. Try substepping.")
success = subStepAnalyze(dT / 2.0, subStep + 1)
curStep = int((curTime - 20000) / dT + 1)
remStep = int(numStep - curStep)
print(f"Current step: {curStep}, Remaining steps: {remStep}")

endT = time.time()
print(f"Loading analysis execution time: {endT - startT} seconds.")

op.wipe()

#%%PostProcessing
import pandas as pd
df_stress = pd.read_csv('stress.txt', sep=" ", header=None)
df_strain = pd.read_csv('strain.txt', sep=" ", header=None)
#df_PP = pd.read_csv('CycPP.txt', sep=" ", header=None)

#stress yy = stress xx
Stress_yy = df_stress.iloc[:, 2].to_numpy()*(-1.0) #compression is positive
Principle_Stress = df_stress.iloc[:, 3].to_numpy()*-1.0
Principle_Strain = df_strain.iloc[:, 3].to_numpy()*-1.0

fig = plt.figure(figsize=(10,18))
ax0 = fig.add_subplot(211)
ax0.plot(Principle_Strain, Principle_Stress, label='stress-strain', linewidth=1.2, color='red')
ax0.set_xlabel("Principle strain,%", fontsize=24)
ax0.set_ylabel("Principle stress, kPa", fontsize=24)
ax0.legend(fontsize=24)

ax1 = fig.add_subplot(212)
ax1.plot(Stress_yy, Principle_Stress, label='Stress Path', linewidth=1.2, color='red')
ax1.set_xlabel("Stress yy = xx, kPa", fontsize=24)
ax1.set_ylabel("Principle stress, kPa", fontsize=24)
ax1.legend(fontsize=24)
plt.show()
plt.close()