How to freeze the particles in the imported dense granular packing

Hi all,

This question follows my previous one [1].
After successfully imported a dense packing with a specific range of size, I would like to further isotropically compress the cropped dense packing with additional two boundary walls.
When the isotropic compression is running, I find that all the particles also start to move around in any direction. I don't want this to happen since I want to maintain the initial thickness of the packing and all the particle position. I can imagine the reason behind this is that there are forces stored at the grain contacts whenever it is compacted.
I've tried to use "blockDOFs" to fix the position of all the particles but in this way, the desired state of isotropic compression cannot be reached.
So my question is how do I freeze the position of all the particles without removing their original energy states and still reach the desired state of isotropic compression?

Thanks in advance!

[1] https://answers.launchpad.net/yade/+question/694288

Below are my two MWEs for packing generation and imported dense packing.

####################### Packing generation MWE ###########################

from yade import pack,plot,export
sp=pack.SpherePack()
O.periodic=True
a=5
b=5
c=5
RADIUS_dimension=0.0250
OUT='mypacking'
length=a*(2*RADIUS_dimension)
height=b*(2*RADIUS_dimension)
width=c*(2*RADIUS_dimension)
thickness=RADIUS_dimension

### Particle size distribution
psdSizes=[.0232,.0254,.0276]
psdCumm=[0,0.5,1]
### friction angles
spFRIC=0.5
### boundary conditions
PI=1.e5
#####Mair and Hazzard (2007)'s microproperties
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=spFRIC,normalCohesion=-1,shearCohesion=-1,label='sphereMat'))

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)
sp.makeCloud((0,0,0),(length,3*height,width), psdSizes = psdSizes, psdCumm = psdCumm, periodic=True, seed =1)
O.bodies.append([utils.sphere(s[0],s[1],color=(1,1,0),material='sphereMat') for s in sp])

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI,-PI,-PI),globUpdate=1,maxStrainRate=(10,10,10),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

def triaxDone():
 triax.dead=True
 O.pause()

O.run(10000000,1)
export.text(OUT+".isoCompression")

#################### Isotropic compression MWE #######################

from yade import pack,plot,export,ymport

sp=pack.SpherePack()

O.periodic=True

a=2.5
b=2.5
c=2.5
RADIUS_dimension=0.0250
OUT='mypacking'
length=a*(2*RADIUS_dimension)
height=b*(2*RADIUS_dimension)
width=c*(2*RADIUS_dimension)
thickness=RADIUS_dimension

# Old packing size for spheres generation
a_old=5
b_old=5
c_old=5
length_old=a_old*(2*RADIUS_dimension)
height_old=b_old*(2*RADIUS_dimension)
width_old=c_old*(2*RADIUS_dimension)
thickness_old=RADIUS_dimension

### Packing ratio
packingRatio = a / a_old

### Particle size distribution
psdSizes=[.0232,.0254,.0276]
psdCumm=[0,0.5,1]

# friction angles
wallFRIC=90
spFRIC=0.5

# boundary conditions
PI=1.e5

### Material properties
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=radians(wallFRIC),normalCohesion=-1,shearCohesion=-1,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=spFRIC,normalCohesion=-1,shearCohesion=-1,label='sphereMat'))

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

upBox = utils.box(center=(length/2,2*height+1.5*thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/5.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(length/2,height-1.5*thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/5.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')

O.bodies.append([upBox,lowBox])

### Import isotropic-compressed packing and boundary
packing = ymport.text("mypacking.isoCompression",color=(1,1,0),material='sphereMat')

### Filter only desired packing
def sphereWanted(b):
    x,y,z = b.state.pos
    return x >= 0 and x <= packingRatio*length_old and y >= packingRatio*height_old and y <= 2*packingRatio*height_old and z >= 0 and z <= packingRatio*width_old

packing = [b for b in packing if sphereWanted(b)]

O.bodies.append(packing)

### Fix particle position
#for b in O.bodies:
# if isinstance(b.shape,Sphere):
# b.state.blockedDOFs='xyzXYZ'
# b.state.dynamics = False
# b.state.vel = Vector3(0,0,0)

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

def triaxDone():
 triax.dead=True
 O.pause()

O.run(100000000,1)

print ('Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
print ('Normal stress (contacts) = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1])

#################### END ###############################

Cheers,
Chien-Cheng