Fixed some particles in a periodic cell

Dear all,

I am trying to fix some particles with blockedDOFs, while others could move freely under external force. It works fine in an aperiodic scene, but failed in a periodic scene.

Following is a minimum code to present the situation. I print the relative position of two particles, which are supposed to be fixed. Also, the velocity of them is indeed not 0, although I set them to 0 in the script.

Thanks a lot for any help.
Huanran

####################################################
from yade import pack

coarseNum = 20 # coarse Number
fineNum = 10 # fine number
coarseSize = 2.5e-3
fineSize = 0.4e-3
cellLength = 0.05 # initial cell size
mn,mx=Vector3(0,0,0),Vector3(cellLength, cellLength, cellLength)
sigma0 = -5e3 # 5kPa confining pressure
F = 2350 * 4./3. * pi * fineSize**3 * 10 #force
F *= Vector3(0,-1,0)

O.materials.append(FrictMat(young=8.e6, poisson=0.8, frictionAngle=radians(26.6), density=2350, label='mat1'))

sp=pack.SpherePack()
sp.makeCloud(minCorner=mn, maxCorner=mx, rMean=coarseSize, rRelFuzz=0., num=coarseNum, periodic=True, seed=1)
sp.makeCloud(minCorner=mn, maxCorner=mx, rMean=fineSize, rRelFuzz=0., num=fineNum, periodic=True, seed=1)
sp.toSimulation(material=O.materials['mat1'],color=(0,0,1))
for p in O.bodies:
   if p.shape.radius == fineSize:
      p.shape.color = (1,0,0)

Gl1_Sphere.stripes=1
yade.qt.Controller(), yade.qt.View()

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   GlobalStiffnessTimeStepper(active=1, timeStepUpdateInterval=100, timestepSafetyCoefficient=0.8),
   PeriTriaxController(
      dynCell=True,
      goal=(sigma0,sigma0,sigma0),
      stressMask=7,
      absStressTol=1.,
      relStressTol=1.e-4,
      maxUnbalanced=1.e-4,
      maxStrainRate=(1.,1.,1.),
      doneHook='applyForce()',
      label='triax'
   ),
   NewtonIntegrator(damping=.3,label='newton')
]

def applyForce():
   triax.dead = True
   O.cell.velGrad = Matrix3(0,0,0, 0,0,0, 0,0,0)
   for p in O.bodies:
      if p.shape.radius > fineSize:
         p.state.blockedDOFs = 'xyzXYZ'
         p.state.vel = Vector3(0,0,0)
         p.state.angVel = Vector3(0,0,0)
      if p.shape.radius == fineSize:
         O.forces.setPermF(p.id, F)
   O.engines+=[PyRunner(command='print(O.bodies[0].state.pos-O.bodies[1].state.pos)',iterPeriod=100,nDo=20)]

print("computing, be patient...")

O.run()