Use 9_4_quadUP element but python stops working

[xzzgameboy]

I am practicing a 1-D consolidation with 9_4_quadUP elements.

The original tcl code can be found at https://opensees.berkeley.edu/wiki/inde ... solidation

The problem I have is that

op.element('9_4_QuadUP', eleTag, nI, nJ, nK, nL, nM, nN, nP, nQ, nR, thick, matTag, bulk, H2ODensity, hPerm, vPerm, wgtX, wgtY)

this line always leads to the restart of python kernel for some reason I do not understand.

The code runs until the line with element command.



Python version: 3.9.12


Any help is greatly appreciated.

Regards,

Steve

Python Code:

Code: Select all

[
from datetime import datetime
import openseespy.opensees as op
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["font.family"] = "Times New Roman"
import sys
#%%
# specify desired drainage type; 1 --> single drainage, 2 --> double drainage
drainage = 1

if drainage == 1:
    print('Single Drainage')
elif drainage == 2:
    print('Double Drainage')
else:
    print('Wrong drainage condition!')
    sys.exit(0)

#%%
#create model
#Remove the existing model, important!!! 
op.wipe()

#%%
#=========================================
#1. Define Mesh Geometry
#=========================================

# number of elements in horizontal direction
numEleX = 1

# size of elements in horizontal direction (m)
sizeEleX = 1.0

# number of nodes in horizontal direction
numNodeX = numEleX* 2 + 1

# number of elements in vertical direction
numEleY = 11  # numEleY should be an odd number
# size of elements in vertical direction (m)
sizeEleY = 1.0
# number of nodes in vertical direction 
numNodeY = 2 * numEleY+1

# total number of nodes
totalNumNode = 3*numNodeY
# total number of elements
totalNumEle = numEleY*numEleX

node_list = []
#-------------------------------------------------------------------------------------------
#  2. DEFINE CORNER NODES
#-------------------------------------------------------------------------------------------
op.model('basic', '-ndm', 2, '-ndf', 3) #corner nodes record PP
fig = plt.figure(figsize=(5,10))
ax0 = fig.add_subplot(211)
ax0.set_xlabel("x, m", fontsize=16)
ax0.set_ylabel("y, m", fontsize=16)
ax0.set_aspect('equal')
ax0.set_xlim(-1.0, 2.0)
ax0.set_ylim(-1.0,numEleY+1)
for i in range(numEleY):
    if (i % 2) == 0: 
        nodeNum0 = 6*i
        nodeNum1 = 6*i+2
        nodeNum2 = 6*i+6
        nodeNum3 = 6*i+8
        node_list.append(nodeNum0)
        node_list.append(nodeNum1)
        node_list.append(nodeNum2)
        node_list.append(nodeNum3)
        Xcoord0 = 0.0
        Xcoord1 = 1.0
        Xcoord2 = 0.0
        Xcoord3 = 1.0
        Ycoord0 = i*1.0
        Ycoord1 = i*1.0
        Ycoord2 = i*1.0 + 1.0
        Ycoord3 = i*1.0 + 1.0
        #set corner nodes
        op.node(nodeNum0, Xcoord0, Ycoord0)
        op.node(nodeNum1, Xcoord1, Ycoord1)
        op.node(nodeNum2, Xcoord2, Ycoord2)
        op.node(nodeNum3, Xcoord3, Ycoord3)
        #plot corner nodes
        ax0.plot(Xcoord0, Ycoord0, 'ro', label = str(nodeNum0))
        ax0.plot(Xcoord1, Ycoord1, 'ro', label = str(nodeNum1))
        ax0.plot(Xcoord2, Ycoord2, 'ro', label = str(nodeNum2))
        ax0.plot(Xcoord3, Ycoord3, 'ro', label = str(nodeNum3))

print('Finished creating corner soil nodes...')
nodeNum_surface1 = nodeNum3  #corner node on the surface
nodeNum_surface2 = nodeNum2
print('Number of nodes are:',nodeNum3+1)
#-------------------------------------------------------------------------------------------
#  3. DEFINE BOUNDARY CONDITIONS FOR CORNER NODES
#-------------------------------------------------------------------------------------------

# boundary conditions for base nodes depending upon selected drainage type


if drainage == 1:       #single drainage, base fixed
    op.fix(0, 1, 1, 0)
    op.fix(2, 1, 1, 0)
elif drainage == 2:       #double drainages, base deformation fixed, drainage free
    op.fix(0, 1, 1, 1)
    op.fix(2, 1, 1, 1)    

#change the surface node to free-drainage no displacement in both directions
op.fix(nodeNum_surface1, 1 , 1, 0)
op.fix(nodeNum_surface2, 1 , 1, 0)

#print('Nodes with fixed X and Y displacement are:')
for i in range(numEleY): #first fix all other element in the x direction
    if (i % 2) == 0: 
        nodeNum0 = 6*i
        nodeNum1 = 6*i+2
        nodeNum2 = 6*i+6
        nodeNum3 = 6*i+8
        if i != 0:
            op.fix(nodeNum0, 1, 0, 0)
            op.fix(nodeNum1, 1, 0, 0)
            #print(nodeNum0)
            #print(nodeNum1)
            
        if i != numEleY-1:
            op.fix(nodeNum2, 1, 0, 0)
            op.fix(nodeNum3, 1, 0, 0)
            #print(nodeNum2)
            #print(nodeNum3)
        

print('Finished creating boundary conditions for corner nodes...')

#-------------------------------------------------------------------------------------------
#  4. DEFINE INTERIOR NODES
#-------------------------------------------------------------------------------------------
op.model('basic', '-ndm', 2, '-ndf', 2) #interior nodes do not record PP, 2 dof
for i in range(numEleY):
    nodeNum0 = 6*i+1
    nodeNum1 = 6*i+3
    nodeNum2 = 6*i+4
    nodeNum3 = 6*i+5
    Xcoord0 = 0.5
    Xcoord1 = 0.0
    Xcoord2 = 0.5
    Xcoord3 = 1.0
    Ycoord0 = i*1.0
    Ycoord1 = i+0.5
    Ycoord2 = i+0.5
    Ycoord3 = i+0.5
    node_list.append(nodeNum0)
    node_list.append(nodeNum1)
    node_list.append(nodeNum2)
    node_list.append(nodeNum3)
    #set interior nodes
    op.node(nodeNum0, Xcoord0, Ycoord0)
    op.node(nodeNum1, Xcoord1, Ycoord1)
    op.node(nodeNum2, Xcoord2, Ycoord2)
    op.node(nodeNum3, Xcoord3, Ycoord3)
    ax0.plot(Xcoord0, Ycoord0, 'bo', label = str(nodeNum0))
    ax0.plot(Xcoord1, Ycoord1, 'bo', label = str(nodeNum1))
    ax0.plot(Xcoord2, Ycoord2, 'bo', label = str(nodeNum2))
    ax0.plot(Xcoord3, Ycoord3, 'bo', label = str(nodeNum3))

nodeNum_surface3 = 6*i+7
Xcoord = 0.5
Ycoord = i+1.0
op.node(nodeNum_surface3, Xcoord, Ycoord)
ax0.plot(Xcoord, Ycoord, 'bo', label = str(nodeNum0))
print('Finished creating all soil nodes...')

#-------------------------------------------------------------------------------------------
#  5. DEFINE BOUNDARY CONDITIONS FOR INTERIOR NODES
#-------------------------------------------------------------------------------------------

op.fix(1, 1, 1)
for i in range(numEleY):
    nodeNum0 = 6*i+1
    nodeNum1 = 6*i+3
    nodeNum2 = 6*i+4
    nodeNum3 = 6*i+5
    if i==0:
        op.fix(nodeNum1, 1, 0)
        op.fix(nodeNum2, 1, 0)
        op.fix(nodeNum3, 1, 0)
        #print(nodeNum1,nodeNum2,nodeNum3)
    else:
        op.fix(nodeNum0, 1, 0)
        op.fix(nodeNum1, 1, 0)
        op.fix(nodeNum2, 1, 0)
        op.fix(nodeNum3, 1, 0)
        #print(nodeNum0,nodeNum1,nodeNum2,nodeNum3)
op.fix(nodeNum_surface3, 1, 0)

#make surface node equal displacement at  degrees 1 & 2
op.equalDOF(nodeNum_surface1,nodeNum_surface2,1,2)
op.equalDOF(nodeNum_surface1,nodeNum_surface3,1,2)
print('Finished creating all boundary conditions and equalDOF...')

#-------------------------------------------------------------------------------------------
#  6. DEFINE MATERIAL PROPERTIES OF SOIL

#-------------------------------------------------------------------------------------------

# Constutive Model: PressureIndependMultiYield material


# saturated mass density (Mg/m3)
satDensity = 2.3
# fluid mass density (Mg/m3)
H2ODensity = 1.0

# shear modulus (kPa)
shear = 2.5e4
# bulk modulus (kPa)
bulk = 6.2e5
# undrained bulk modulus (kPa)
uBulk = 2.2e5

# vertical permeability, input value in m/s, units are modified to suit constitutive model
mPermV = 5.0e-5
vPerm = mPermV/9.81/H2ODensity

# horizontal permeability, input value in m/s, units are modified to suit constitutive model
mPermH = 5.0e-5
hPerm = mPermH/9.81/H2ODensity

# cohesion (kPa)
cohesion = 45.0

# friction angle (degrees)
phi = 0.0

# peak shear strain 
gammaPeak = 0.1

# reference pressure (kPa)
refP = 80.0

# pressure dependency coefficient
pressCoeff = 0.0

# number of yield surfaces
numSurf = 22

#-------------------------------------------------------------------------------------------
#  7. DEFINE SOIL MATERIALS
#-------------------------------------------------------------------------------------------
'''
Soil Model: PressureIndependMultiYield 
nDMaterial('PressureIndependMultiYield', matTag, nd, rho, refShearModul, refBulkModul, 
           cohesi, peakShearStra, frictionAng=0., refPress=100., pressDependCoe=0., 
           noYieldSurf=20, *yieldSurf)
'''

matTag = 1

nd = 2 #Number of dimensions, 2 for plane-strain, and 3 for 3D analysis.
rho = satDensity # Saturated Soil Mass Density

op.nDMaterial('PressureIndependMultiYield', matTag, 2, satDensity, shear, bulk, 
              cohesion, gammaPeak, phi, refP, 0., numSurf)
print('Finished creating all soil materials....')

#-------------------------------------------------------------------------------------------
#  8. DEFINE SOIL ELEMENTS
#-------------------------------------------------------------------------------------------
#element('9_4_QuadUP', eleTag, *eleNodes, thick, matTag, bulk, fmass, hPerm, vPerm, <b1=0, b2=0>)


thick = 1.0
wgtX = 0.0
wgtY = -9.81
i = 0
n0 = 6*i
n1 = n0+1
n2 = n0+2
n3 = n0+3
n4 = n0+4
n5 = n0+5
n6 = n0+6
n7 = n0+7
n8 = n0+8

op.element('9_4_QuadUP', i, n0, n1, n2, n3, n4, n5, n6, n7, n8, thick, matTag, bulk, H2ODensity, hPerm, vPerm, wgtX, wgtY)

[xzzgameboy]

sorry that the original code I posted has some errors but I have update the code.

If there is wrong with the code, error or warning should be given by python but the only thing I can see is the warning that kernel dead.

I tried with Spyder and Jupyter Notebook. Both Macos and Windows 11. All same issue.

Any help is greatly appreciated!

Thanks in advance.

Regards

Steve

[mhscott]

What version of OpenSeesPy are you using? version 3.4.0.2 is the newest

[xzzgameboy]

What version of OpenSeesPy are you using? version 3.4.0.2 is the newest

Yes 3.4.0.2, the latest version of openseespy

[mhscott]

You have an issue with the nodal dofs of the interior/exterior nodes.
This error comes out of the element setDomain method:
"FATAL ERROR NineFourNodeQuadUP, has wrong number of DOFs at its nodes 0"
So, this is not an OpenSeesPy issue. Make sure you converted the Tcl script to Python correctly.

[mhscott]

With a little more useful error messages, it seems the issue is with node 1
"FATAL ERROR NineFourNodeQuadUP tag=0 has wrong number of DOFs at node 1"

[xzzgameboy]

With a little more useful error messages, it seems the issue is with node 1
"FATAL ERROR NineFourNodeQuadUP tag=0 has wrong number of DOFs at node 1"

Thank you for all the help Prof. Scott! I will try to correct the error and put the correct code here.

[xzzgameboy]

Hi all,

Just in case someone is also interested in this topic.

The reason my code didn't work was the I set the wrong boundary condition on two nodes (dof 3, pore pressure).

Code: Select all

[/
# -*- coding: utf-8 -*-
"""
Created on Sat Jul 23 11:13:25 2022

@author: Zhongze Xu, the University of Texas at Austin

1D consolidation using fully coupled, u-p formulation.
Cohesive soil. Option of single or double drainage. 
Basic units are kN, m, and sec (unless specified)


Details of the model can be found:
https://opensees.berkeley.edu/wiki/index.php/One-dimensional_Consolidation

Original tcl script was written by:
    Christopher McGann and Pedro Arduino, University of Washington
"""
#from IPython import get_ipython;   
#get_ipython().run_line_magic('reset', '-sf')

from datetime import datetime
import openseespy.opensees as op
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["font.family"] = "Times New Roman"
import sys
from matplotlib import cm
import pandas as pd
#%%
# specify desired drainage type; 1 --> single drainage, 2 --> double drainage
drainage = 1

if drainage == 1:
    print('Single Drainage')
elif drainage == 2:
    print('Double Drainage')
else:
    print('Wrong drainage condition!')
    sys.exit(0)

#%%
#create model
#Remove the existing model, important!!! 
op.wipe()

#%%
#=========================================
#1. Define Mesh Geometry
#=========================================

# number of elements in horizontal direction
numEleX = 1

# size of elements in horizontal direction (m)
sizeEleX = 1.0

# number of nodes in horizontal direction
numNodeX = numEleX* 2 + 1

# number of elements in vertical direction
numEleY = 10
# size of elements in vertical direction (m)
sizeEleY = 1.0
# number of nodes in vertical direction 
numNodeY = 2 * numEleY+1

# total number of nodes
totalNumNode = 3*numNodeY
# total number of elements
totalNumEle = numEleY*numEleX

node_list = []
#%%
#-------------------------------------------------------------------------------------------
#  2. DEFINE CORNER NODES
#-------------------------------------------------------------------------------------------
op.model('basic', '-ndm', 2, '-ndf', 3) #corner nodes record PP
fig = plt.figure(figsize=(10,10))
ax0 = fig.add_subplot(121)
ax0.set_xlabel("x, m", fontsize=16)
ax0.set_ylabel("y, m", fontsize=16)
ax0.set_aspect('equal')
ax0.set_xlim(-1.0, 2.0)
ax0.set_ylim(-1.0,numEleY+1)
for i in range(0,numEleX+2,2):
    #print(i)
    for j in range(0,2*numEleY+2,2):
        #print(j)
        nodeNum = i+1 + (j)*numNodeX
        node_list.append(nodeNum)
        Xcoord = (i)*sizeEleX/2
        Ycoord = (j)*sizeEleX/2
        #set corner nodes
        op.node(nodeNum, Xcoord, Ycoord)

        #plot corner nodes
        ax0.plot(Xcoord, Ycoord, 'ro', label = str(nodeNum))


print('Finished creating corner soil nodes...')

#%%
#-------------------------------------------------------------------------------------------
#  3. DEFINE BOUNDARY CONDITIONS FOR CORNER NODES
#-------------------------------------------------------------------------------------------

# boundary conditions for base nodes depending upon selected drainage type


if drainage == 1:       #single drainage, base deformation fixed, pore pressure free = no drainage
    op.fix(1, 1, 1, 0)
    op.fix(3, 1, 1, 0)
elif drainage == 2:       #double drainages, base deformation fixed, pore pressure fixed, free drainage 
    op.fix(1, 1, 1, 1)
    op.fix(3, 1, 1, 1)    

#change the surface node to free-drainage no displacement in both directions
op.fix(totalNumNode, 1 , 0, 1)
op.fix(totalNumNode-2, 1 , 0, 1)

#print('Nodes with fixed X and Y displacement are:')

for i in range(numEleY): #first fix all other element in the x direction
    if i != 0:
        nodeNum3 = 6*(i+1)-5
        nodeNum4 = 6*(i+1)-3
        op.fix(nodeNum3, 1, 0, 0)
        op.fix(nodeNum4, 1, 0, 0)  
        print(nodeNum3)
        print(nodeNum4)
        

print('Finished creating boundary conditions for corner nodes...')
#%%
#-------------------------------------------------------------------------------------------
#  4. DEFINE INTERIOR NODES
#-------------------------------------------------------------------------------------------
op.model('basic', '-ndm', 2, '-ndf', 2) #interior nodes do not record PP, 2 dof
for i in range(numEleY):
    nodeNum0 = 6*(i+1)-4
    nodeNum1 = 6*(i+1)-2
    nodeNum2 = 6*(i+1)-1
    nodeNum3 = 6*(i+1)
    Xcoord0 = 0.5
    Xcoord1 = 0.0
    Xcoord2 = 0.5
    Xcoord3 = 1.0
    Ycoord0 = i*1.0
    Ycoord1 = i+0.5
    Ycoord2 = i+0.5
    Ycoord3 = i+0.5
    node_list.append(nodeNum0)
    node_list.append(nodeNum1)
    node_list.append(nodeNum2)
    node_list.append(nodeNum3)
    print(nodeNum0, nodeNum1, nodeNum2, nodeNum3)
    #set interior nodes
    op.node(nodeNum0, Xcoord0, Ycoord0)
    op.node(nodeNum1, Xcoord1, Ycoord1)
    op.node(nodeNum2, Xcoord2, Ycoord2)
    op.node(nodeNum3, Xcoord3, Ycoord3)
    ax0.plot(Xcoord0, Ycoord0, 'bo', label = str(nodeNum0))
    ax0.plot(Xcoord1, Ycoord1, 'bo', label = str(nodeNum1))
    ax0.plot(Xcoord2, Ycoord2, 'bo', label = str(nodeNum2))
    ax0.plot(Xcoord3, Ycoord3, 'bo', label = str(nodeNum3))

nodeNum_surface3 = 6*(i+1)+2
Xcoord = 0.5
Ycoord = i+1.0
op.node(nodeNum_surface3, Xcoord, Ycoord)
ax0.plot(Xcoord, Ycoord, 'bo', label = str(nodeNum0))
print('Finished creating all soil nodes...')
#%%
#-------------------------------------------------------------------------------------------
#  5. DEFINE BOUNDARY CONDITIONS FOR INTERIOR NODES
#-------------------------------------------------------------------------------------------

op.fix(2, 1, 1)
for i in range(numEleY):
    nodeNum0 = 6*(i+1)-4
    nodeNum1 = 6*(i+1)-2
    nodeNum2 = 6*(i+1)-1
    nodeNum3 = 6*(i+1)
    if i==0:
        op.fix(nodeNum1, 1, 0)
        op.fix(nodeNum2, 1, 0)
        op.fix(nodeNum3, 1, 0)
        #print(nodeNum1,nodeNum2,nodeNum3)
    else:
        op.fix(nodeNum0, 1, 0)
        op.fix(nodeNum1, 1, 0)
        op.fix(nodeNum2, 1, 0)
        op.fix(nodeNum3, 1, 0)
        #print(nodeNum0,nodeNum1,nodeNum2,nodeNum3)
op.fix(nodeNum_surface3, 1, 0)

nodeNum_surface1 = 6*numEleY+3
nodeNum_surface2 = 6*numEleY+1

#make surface node equal displacement at  degrees 1 & 2
op.equalDOF(nodeNum_surface1,nodeNum_surface2,1,2)
op.equalDOF(nodeNum_surface1,nodeNum_surface3,1,2)
print('Finished creating all boundary conditions and equalDOF...')
plt.show()
plt.close()
#-------------------------------------------------------------------------------------------
#  6. DEFINE MATERIAL PROPERTIES OF SOIL

#-------------------------------------------------------------------------------------------

# Constutive Model: PressureIndependMultiYield material


# saturated mass density (Mg/m3)
satDensity = 2.3
# fluid mass density (Mg/m3)
H2ODensity = 1.0

# shear modulus (kPa)
shear = 2.5e4
# bulk modulus (kPa)
bulk = 6.2e5
# undrained bulk modulus (kPa)
uBulk = 2.2e5

# vertical permeability, input value in m/s, units are modified to suit constitutive model
mPermV = 5.0e-5
vPerm = mPermV/9.81/H2ODensity

# horizontal permeability, input value in m/s, units are modified to suit constitutive model
mPermH = 5.0e-5
hPerm = mPermH/9.81/H2ODensity

# cohesion (kPa)
cohesion = 45.0

# friction angle (degrees)
phi = 0.0

# peak shear strain 
gammaPeak = 0.1

# reference pressure (kPa)
refP = 80.0

# pressure dependency coefficient
pressCoeff = 0.0

# number of yield surfaces
numSurf = 22

#-------------------------------------------------------------------------------------------
#  7. DEFINE SOIL MATERIALS
#-------------------------------------------------------------------------------------------
'''
Soil Model: PressureIndependMultiYield 
nDMaterial('PressureIndependMultiYield', matTag, nd, rho, refShearModul, refBulkModul, 
           cohesi, peakShearStra, frictionAng=0., refPress=100., pressDependCoe=0., 
           noYieldSurf=20, *yieldSurf)
'''

matTag = 1

nd = 2 #Number of dimensions, 2 for plane-strain, and 3 for 3D analysis.
rho = satDensity # Saturated Soil Mass Density

op.nDMaterial('PressureIndependMultiYield', matTag, 2, satDensity, shear, bulk, 
              cohesion, gammaPeak, phi, refP, 0., numSurf)
print('Finished creating all soil materials....')

#-------------------------------------------------------------------------------------------
#  8. DEFINE SOIL ELEMENTS
#-------------------------------------------------------------------------------------------
#element('9_4_QuadUP', eleTag, *eleNodes, thick, matTag, bulk, fmass, hPerm, vPerm, <b1=0, b2=0>)


thick = 1.0
wgtX = 0.0
wgtY = -9.81
'''
i = 0
n0 = 6*i
n1 = n0+1
n2 = n0+2
n3 = n0+3
n4 = n0+4
n5 = n0+5
n6 = n0+6
n7 = n0+7
n8 = n0+8
#op.element('9_4_QuadUP', 1, n0, n1, n2, n3, n4, n5, n6, n7, n8, thick, matTag, bulk, H2ODensity, hPerm, vPerm, wgtX, wgtY)
'''


for i in range(numEleY):
    eleTag = i+1
    n1 = 6*(i+1)-5
    n2 = n1+1
    n3 = n1+2
    n4 = n1+3
    n5 = n1+4
    n6 = n1+5
    n7 = n1+6
    n8 = n1+7
    n9 = n1+8
    print('Node:', n1, n3, n9, n7, n2, n6, n8, n4, n5) 
    #Need to check the order of nodes input to the element commend\
    #https://opensees.berkeley.edu/wiki/index.php/Four_Node_Quad_u-p_Element
    op.element('9_4_QuadUP', eleTag, n1, n3, n9, n7, n2, n6, n8, n4, n5, thick, 
               matTag, bulk, H2ODensity, hPerm, vPerm, wgtX, wgtY)



# update material stage to zero for elastic behavior
op.updateMaterialStage('-material', matTag, '-stage', 0)

print('Finished creating all soil elements...')
#%%
#-------------------------------------------------------------------------------------------
#  8. APPLY GRAVITY LOADING
#-------------------------------------------------------------------------------------------

gamma = 0.5
beta = 0.25

#Analysis officially starts here
op.constraints('Transformation')
op.test('NormDispIncr', 1.0e-5, 30, 1)
op.algorithm('Newton')
op.numberer('RCM')
op.system('ProfileSPD')
op.integrator('Newmark', gamma, beta)
op.analysis('Transient')

op.analyze(10, 500)

print('Finished with elastic gravity analysis...')

#%%
#-------------------------------------------------------------------------------------------
#  9. UPDATE SOIL MATERIAL STAGE
#-------------------------------------------------------------------------------------------

# update material stage to consider elastoplastic behavior

op.updateMaterialStage('-material', matTag, '-stage', 1)
op.analyze(40, 500)

print('Finished with plastic gravity analysis...')

#%%
#-------------------------------------------------------------------------------------------
#  10. CREATE RECORDERS
#-------------------------------------------------------------------------------------------
op.recorder('Node','-file', 'displacement.out','-time', '-nodeRange',1,totalNumNode,'-dof', 1, 2, 'disp')
op.recorder('Node','-file', 'porepressure.out','-time', '-nodeRange',1,totalNumNode,'-dof', 3, 'vel')

#record stress and strain at the 9th Gaussian point, center of the element
op.recorder('Element','-file','stress.out','-time','-eleRange', 1, numEleY,'material','9','stress')
op.recorder('Element','-file','strain.out','-time','-eleRange', 1, numEleY,'material','9','strain')

print('Finished creating all recorders...')

#%%
#-------------------------------------------------------------------------------------------
#  11. CONSOLIDATION ANALYSIS  
#-------------------------------------------------------------------------------------------

setTime = 0.0

op.wipeAnalysis()

overburden  = -1000.0 #set vertical load 

op.timeSeries('Constant', 1) #tsTag = 1
op.pattern('Plain', 1, 1)
op.load(totalNumNode-1, 0.0,overburden) #put overburden on y direction

print('Consolidation loading created...')

op.constraints('Penalty', 1.0e16, 1.0e16)
op.test('NormDispIncr', 1.0e-4, 30, 1)
op.algorithm('Newton')
op.numberer('RCM')
op.system('ProfileSPD')
op.integrator('Newmark', gamma, beta)
op.rayleigh(0.5, 0.5, 0.0, 0.0)
op.analysis('Transient')

if drainage == 1:
    op.analyze(800, 0.25)
else:
    op.analyze(2000, 0.25)

print('Finished with consolidation analysis...')

op.wipe()

#%%
#-------------------------------------------------------------------------------------------
#  12. Post-Process
#-------------------------------------------------------------------------------------------
df = pd.read_csv('porepressure.out', delimiter=' ', header=None)

df_ID = df.iloc[:,0]
step_ID = df_ID.to_numpy()
step_ID = step_ID-step_ID[0]

df =  df.iloc[: , 1:]

nStep = df.shape[0]
nNode = df.shape[1]
plot_times = nStep
DOC = step_ID/step_ID[-1]  #something like degree of consolidation for plot only!!!!
depth_left = np.linspace(-numEleY, 0, numEleY+1)

fig= plt.figure(figsize=(6,6))
ax0 = fig.add_subplot(111)
ax0.set_xlabel('Excess Pore Pressure, kPa', fontsize=20,  fontweight='bold')
ax0.set_ylabel('Depth, m', fontsize=20, fontweight='bold')
ax0.set_ylim(-numEleY-0.5,0.5)
ax0.tick_params(direction='in', labelsize=13)
ax0.tick_params(axis='both', which='both',direction='in',labelsize=14, width=1.5, length=3)


for i in range(plot_times):
    plot_arr = df.iloc[i,:].to_numpy()
    PP_left = plot_arr[0:nNode:6]
    DOC_arr = DOC[i]*np.ones(numEleY+1)
    if i > 5:
        im = ax0.plot(PP_left,depth_left, lw=3, c=cm.turbo(DOC[i]))

if drainage == 1:
    ax0.text(0.9, 0.9,'Single Drainage', fontsize = 18,horizontalalignment='right',
             verticalalignment='top',transform = ax0.transAxes)
else:
    ax0.text(0.9, 0.9,'Double Drainage', fontsize = 18,horizontalalignment='right',
             verticalalignment='top',transform = ax0.transAxes)
plt.show()
plt.close()

[mhscott]

Thanks for sharing the solution!