Some Problems about Direct Shear Test

Hello, everyone. I'm trying to simulate the direct shear process with the following program[1]. Before running the script, I generated particles in another script. But when I run the following script, the system can't find the particles. But when I input the commands of the script sentence by sentence in the Yade program (before importing the particles), I can import the particles normally. Is there someone can tell me where the problem is?
[1]the script:
from yade import ymport, utils, plot
from yade import ymport, utils, plot
from yade import pack,export
from yade import polyhedra_utils
from yade import geom,utils,plot
from yade import qt,ymport,timing
from yade import *

################## parameters definition
DAMPING=0.5
dtCOEFF=0.5
normalSTRESS=5e5
normalVEL=normalSTRESS/1e8 # 0.001 for 100kPa , i dont know the reason
shearVEL=1*normalVEL # try different values to ensure quasi-static conditions
intR=1.263 #1.263 for X1Y1Z1_5k
DENS=3000
YOUNG=50e9
FRICT=18
ALPHA=0.3
TENS=45e5
COH=45e6
iterMAX=400000
###################### Import of the sphere assembly
### material definition
sphereMat = FrictMat()
sphereMat.young=YOUNG
sphereMat.frictionAngle=radians(FRICT)
sphereMat.density=DENS

wallMat = FrictMat()
wallMat.young=YOUNG
wallMat.frictionAngle=radians(0)
wallMat.density=DENS

O.materials.append(sphereMat)
O.materials.append(wallMat)

## copy spheres from the packing into the scene
O.bodies.append(ymport.textExt('/tmp/compressed.txt',format='x_y_z_r',material=sphereMat))
## preprocessing to get dimensions
dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf
## initial surface
S0=X*Z
## spheres factory
R=0
Rmax=0
numSpheres=0
for o in O.bodies:
 if isinstance(o.shape,Sphere):
  o.shape.color=(0.7,0.5,0.3)
  numSpheres+=1
  R+=o.shape.radius
  if o.shape.radius>Rmax:
   Rmax=o.shape.radius
  Rmean=R/numSpheres
###################### creation of shear box
thickness=Y/100
oversizeFactor=1.05
O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,3*(Y/4)+yinf,zsup*3/4),extents=(thickness/2,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
butee=O.bodies[-1]
### add top box
O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,3*(Y/4)+yinf,zsup*3/4),extents=(thickness/2,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
O.bodies.append(utils.box(center=(xinf+X/2,3*(Y/4)+yinf,zinf-thickness+Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness/2),material=wallMat,fixed=True,wire=True))
O.bodies.append(utils.box(center=(xinf+X/2,3*(Y/4)+yinf,zsup+thickness-Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness/2),material=wallMat,fixed=True,wire=True))
### add bottom box
O.bodies.append(utils.box(center=(xsup+thickness-Rmean/10,(Y/4)+yinf,zsup*3/4),extents=(thickness/2,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
piston1=O.bodies[-1]
O.bodies.append(utils.box(center=(xinf-thickness+Rmean/10,(Y/4)+yinf,zsup*3/4),extents=(thickness/2,Y/4,oversizeFactor*Z/2),material=wallMat,fixed=True))
piston2=O.bodies[-1]
O.bodies.append(utils.box(center=(xinf+X/2,(Y/4)+yinf,zinf-thickness+Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness/2),material=wallMat,fixed=True,wire=True))
O.bodies.append(utils.box(center=(xinf+X/2,(Y/4)+yinf,zsup+thickness-Rmean/10),extents=(oversizeFactor*X/2,Y/4,thickness/2),material=wallMat,fixed=True,wire=True))
### loading platens,it is a xz plate?
O.bodies.append(utils.box(center=(xinf+X/2.,yinf-thickness+Rmean/10.,zinf+Z/2.),extents=(oversizeFactor*X/2,thickness/2,oversizeFactor*Z/2),material=wallMat,fixed=True))
bottomPlate=O.bodies[-1]
O.bodies.append(utils.box(center=(xinf+X/2.,ysup+thickness-Rmean/10.,zinf+Z/2.),extents=(oversizeFactor*X/2,thickness/2,oversizeFactor*Z/2),material=wallMat,fixed=True))
topPlate=O.bodies[-1]

###################### Engines
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
  ),
 GlobalStiffnessTimeStepper(defaultDt=0.1*PWaveTimeStep(),timestepSafetyCoefficient=dtCOEFF),
 TranslationEngine(ids=[topPlate.id],translationAxis=(0.,-1.,0.),velocity=0.,label='yTranslation'),
 PyRunner(iterPeriod=1,initRun=True,command='servoController()',label='servo'),
 NewtonIntegrator(damping=DAMPING,gravity=(0,0,0),label='damper'),
 PyRunner(iterPeriod=iterMAX/1000,initRun=True,command='dataCollector()'),
]
###################### Engines definition ( servoController() and dataCollector() )
shearing=False
sigmaN=0
tau=0
Fs1=0
Fs2=0
Xdispl=0
px0=0
Ydispl=0
py0=topPlate.state.pos[1]
prevTranslation=0
n=0

def servoController():
 global px0, py0, sigmaN, n, Fn1, Fn2, shearing, butee, piston1, piston2
 Fn1=abs(O.forces.f(topPlate.id)[1])
 Fn2=abs(O.forces.f(bottomPlate.id)[1])
 sigmaN=Fn1/(S0-2*Xdispl*Z) # necessary?
 if shearing==False:
  if yTranslation.velocity<normalVEL:
   yTranslation.velocity+=normalVEL/1000
  if sigmaN>(0.9*normalSTRESS):
   yTranslation.velocity=normalVEL*((normalSTRESS-sigmaN)/normalSTRESS)
 if shearing==False and abs((normalSTRESS-sigmaN)/normalSTRESS)<0.001 and unbalancedForce()<0.01:
  yTranslation.velocity=0
  n+=1
 if n>1000 and abs((sigmaN-normalSTRESS)/normalSTRESS)<0.001:
  print 'stress on joint plane = ', utils.forcesOnPlane((X/2,Y/2,Z/2),(0,1,0))/S0

### add engine and initialisation
O.engines=O.engines[:6]+ [TranslationEngine(ids=[piston1.id,piston2.id],translationAxis=(-1.,0.,0.),velocity=0.)]+O.engines[6:]
px0=piston1.state.pos[0]
py0=topPlate.state.pos[1]
shearing=True
print 'shearing! || iteration=', O.iter
if shearing==True:
 #yTranslation.velocity=0. # if you want a constant normal diaplacement control
 yTranslation.velocity=50*normalVEL*((normalSTRESS-sigmaN)/normalSTRESS) # not good enough because not general (coeff needs to be adjusted)
 if O.engines[6].velocity<shearVEL:
  O.engines[6].velocity+=(shearVEL/1000)

def dataCollector():
 global Xdispl, Ydispl, tau, Fs1, Fs2
 if shearing:
  Fs1=abs(O.forces.f(butee.id)[0])
  Fs2=abs(O.forces.f(piston1.id)[0])
  Xdispl=abs(piston1.state.pos[0]-px0)
 tau=Fs1/(S0-2*Xdispl*Z)
 Ydispl=abs(topPlate.state.pos[1]-py0)
 numberOfCohesiveBonds=0
 for i in O.interactions:
  if i.phys.isCohesive:
   numberOfCohesiveBonds+=1
 yade.plot.addData({'t':O.time,'i':O.iter,'Xdispl':Xdispl,'sigmaN':sigmaN,'Fn1':Fn1,'Fn2':Fn2,'tau':tau,'Fs1':Fs1,'Fs2':Fs2,'Ydispl':Ydispl,'unbF':utils.unbalancedForce(),'Ek':utils.kineticEnergy(),'bonds':numberOfCohesiveBonds})
 plot.saveGnuplot(Output)

# defines window plot
plot.plots={'i':('sigmaN','tau')}

################# graphical intervace
from yade import qt
qt.Controller()
qt.View()
v=qt.Renderer()
v.dispScale=(1,1,1) # displacements are scaled for visualization purpose

################# to manage interaction detection factor during the first timestep
O.step();
################# initializes the interaction detection factor to its default value (new contacts will be created between strictly contacting particles only)
ss2d3dg.interactionDetectionFactor=-1.
aabb.aabbEnlargeFactor=-1.

################# simulation starts here
O.run(iterMAX)