The OpenSees Community

xzzgameboy

mhscott wrote: ↑Wed Feb 08, 2023 3:05 pm Use ops.constraints('Transformation')

Thank you Prof. Scott

xzzgameboy

Hi all,

I am tranferring Opensees Tcl to Openseepy but encountered a problem here.

The code is to model a 9-4 quad u-p element under cyclic load (sine) with PDMY03 model.

The original code can be found https://opensees.berkeley.edu/wiki/inde ... mentdriver
and I am working on the "undrained cyclic" case

the part of tcl i am having issue with is:

Code: Select all

 model basic -ndm 2 -ndf 3
	pattern Plain 4 "Sine 0 [expr $period*1000] $period" {
		load 3 [expr $csr*$vertPress*$xsize*$thickness*0.25] 0 0
		load 4 [expr $csr*$vertPress*$xsize*$thickness*0.25] 0 0
	}
	model basic -ndm 2 -ndf 2
	pattern Plain 5 "Sine 0 [expr $period*1000] $period" {
		load 7 [expr $csr*$vertPress*$xsize*$thickness*0.5] 0
	}
	puts "$Analysis_case"
	while {[expr abs([nodeDisp 7 1]/$ysize)]<$target_shear_strain} {
	# while {[getTime]<0.25}
	analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
	puts "time = [getTime] sec"
	}

so my python script correspondingly looks like this

Code: Select all

op.constraints('Plain')
op.test('NormDispIncr', 1.0e-5, 50, 0)
op.algorithm('KrylovNewton')
op.numberer('RCM')
op.system('ProfileSPD')
op.integrator('Newmark', 5.0/6.0, 4.0/9.0)
op.rayleigh(Mproprayleigh, 0.0, Kinitproprayleigh, 0.0) #modification
op.analysis('VariableTransient')

op.model('basic', '-ndm', 2, '-ndf', 3) 
op.timeSeries('Trig', 3, 0.0, period*1000.0, period, '-factor', 1.0, '-shift', 0.0)
op.pattern('Plain', 3, 3, '-fact', 1.0)
op.sp(4, 1, csr*load_conf)
op.sp(3, 1, csr*load_conf)

op.model('basic', '-ndm', 2, '-ndf', 2)
op.timeSeries('Trig', 4, 0.0, period*1000.0, period, '-factor', 1.0, '-shift', 0.0)
op.pattern('Plain', 4, 4, '-fact', 1.0)
op.sp(7, 1, 2.0*csr*load_conf)


target_shear_strain = 0.03
while abs(op.nodeResponse(3,1,1))/0.1 < target_shear_strain:
    op.analyze(1, deltaT, deltaT/100, deltaT, 20)

I did equalDOF of nodes 3, 4, 7 in DOF 1 and 2, nodes 3 and 4 also equalDOF at degree 3 (fluid pressure)
However, the error i got is:

Code: Select all

WARNING PlainHandler::handle() -  non-homogeneos constraint for node 4 homogeneous constraint assumed
WARNING PlainHandler::handle() -  non-homogeneos constraint for node 3 homogeneous constraint assumed
WARNING PlainHandler::handle() -  non-homogeneos constraint for node 7 homogeneous constraint assumed
WARNING PlainHandler::handle() -  constraint at dof 0 already specified for constrained node in MP_Constraint at node 4
WARNING PlainHandler::handle() -  constraint at dof 0 already specified for constrained node in MP_Constraint at node 7

Any help is highly appreciated,

Thank you.

Regards,

Steve

xzzgameboy

Dear all,

I am trying to use updateMaterials function to update bulk modulus of UCSD soil models in openseespy but it seems there is no corresponding function called updateMaterials like the one in Opensees (https://opensees.berkeley.edu/OpenSees/ ... l/4246.htm). I cannot find it in the openseespy manual (https://openseespydoc.readthedocs.io/en ... index.html).

The tcl script is:

updateMaterials -material $matTag bulkModulus [expr $G*2/3.]

How do I do this in openseespy?

Any help is greatly appreciated.

Regards,

Steve

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()

xzzgameboy

mhscott wrote: ↑Thu Jul 28, 2022 8:39 am 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

mhscott wrote: ↑Wed Jul 27, 2022 8:20 pm 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

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

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

mhscott wrote: ↑Tue May 03, 2022 6:16 am ops.timeSeries('Path',1,'-time',...) #https://openseespydoc.readthedocs.io/en ... athTs.html
ops.pattern('Plain',2,1)
ops.sp(3,...)
ops.sp(4,...)

Thank you very much Prof. Scott!

xzzgameboy

Dear all,

I am transferring the opensees pm4sand example to openseespy.

I am currently stuck on one command:

set ts1 "{Series -time {100 80000 1.0e10} -values {1.0 1.0 1.0} -factor 1}"
eval "pattern Plain 2 $ts1 {
sp 3 2 $vDisp
sp 4 2 $vDisp
}"

Openseespy does not seems to have this Series function or in the timeSeries function, it does not have this ramped load. So what is the equivalent function for this in openseespy. I am quite confused.

Any help or suggestion is greatly appreciated.

Regards,

Steve

The OpenSees Community

Search found 10 matches

Re: Sine Wave load on nodes

Sine Wave load on nodes

cannot find function updateMaterials

Re: Use 9_4_quadUP element but python stops working

Re: Use 9_4_quadUP element but python stops working

Re: Use 9_4_quadUP element but python stops working

Re: Use 9_4_quadUP element but python stops working

Use 9_4_quadUP element but python stops working

Re: Transfer TCL to Openseespy

Transfer TCL to Openseespy