ConnectionFinder problem and Bonding broke after the first frame

According to my previous question (https://answers.launchpad.net/esys-particle/+question/693648), I found out that using the ConnectionFinder to define bonding between particles need to set the MaxDist up to the Minimum radius to make all particles bonded to each other.

> I am very curious that, normally the MaxDist in ConnectionFinder usually had low value (around 0.005 - 1) and can bond all particles together, but in my case it seem to need a large value (around Rmin which is up to 10-40). Why the MaxDist is this large?

Then, after I define all the bond and conducted a 2D biaxial testing, all of bonding just appear properly in the first frame, but in the next frame most of them just disappear (Only a few just left and also broke in next 2-3 step).

First, I think it might be the same with this problem (https://answers.launchpad.net/esys-particle/+question/245754) , so I tried to

> (1.) add linear and rotational damping
> (2.) reduce the timestep size 100 lower than the critical size
> (3.) remove the confining pressure in the beginning

but the bond still broke after the first frame. The only way that I can do to make the bond broke the least is to reduce the Rmin of the particle. However, the number of broken bond still more than half in the second frame without applying any confining pressure.

How do I fix this problem? Do I triggered any bugs?

Thank you for an advance .
Wasin

Here is the script of my biaxial test simulation
####################################

#import the appropriate ESyS-Particle modules:
from esys.lsm import *
from esys.lsm.util import *
from esys.lsm.geometry import *
from WallLoader import WallLoaderRunnable
from ServoWallLoader import ServoWallLoaderRunnable

#Set Parameters
Rmin = 40.0
Rmax = 200.0
Xlength = 4000
Ylength = 3000
YoungMod = 2.00e+9
Poisson = 0.30
Cohes = 20.0e+6
TanPhi = 1.00
MuSta = 0.60
MuDyn = 0.40
Rho = 2500.00
ConfiningP = 10.0e+6

TotalTimestep = 2500000

#instantiate a simulation object:
sim = LsmMpi (numWorkerProcesses = 3, mpiDimList = [3,1,1])

#initialise the neighbour search algorithm:
sim.initNeighbourSearch (
   particleType = "RotSphere",
   gridSpacing = 2.2 * Rmax,
   verletDist = Rmin/5.0
)

#set the number of timesteps and timestep increment:
sim.setNumTimeSteps (TotalTimestep)
sim.setTimeStepSize (1.0e-04)

#enforce two-dimensional computations:
sim.force2dComputations (True)

#specify the spatial domain for the simulation
domain = BoundingBox(Vec3(-0.5 * Xlength,-0.5 * Ylength,0), Vec3( 1.5 * Xlength , 1.5 * Ylength, 0))
sim.setSpatialDomain (domain)

#create a prism of spherical particles:
geoRandomBlock = RandomBoxPacker (
   minRadius = Rmin,
   maxRadius = Rmax,
   cubicPackRadius = 2.2000,
   maxInsertFails = 50000,
   bBox = BoundingBox(
      Vec3(0.0, 0.0, 0.0),
      Vec3(Xlength, Ylength, 0.0)
   ),
   circDimList = [False, False, False],
   tolerance = 1.0000e-05
)
geoRandomBlock.generate()
geoRandomBlock_particles = geoRandomBlock.getSimpleSphereCollection()

#add the particles to the simulation object:
for pp in geoRandomBlock_particles:
    centre = pp.getPosn()
    radius = pp.getRadius()
    Y = centre[1]
    if (Y < Rmin/2): # particle is near the base (tag=2)
        pp.setTag (2)
    elif (Y > (Ylength-Rmax/2)): # particle is near the top (tag=3)
        pp.setTag (3)
    else: # particle is inside the shear cell (tag=1)
        pp.setTag (1)
    sim.createParticle(pp)

#set the density of all particles:
sim.setParticleDensity (
   tag = 1,
   mask = -1,
   Density = Rho
)

#set the density of all particles:
sim.setParticleDensity (
   tag = 2,
   mask = -1,
   Density = Rho
)

#set the density of all particles:
sim.setParticleDensity (
   tag = 3,
   mask = -1,
   Density = Rho
)

#bond particles together with bondTag = 1:
sim.createConnections(
   ConnectionFinder(
      maxDist = Rmin,
      bondTag = 1,
      pList = geoRandomBlock_particles
   )
)

#create a wall at the bottom of the model:
sim.createWall (
   name = "bottom_wall",
   posn = Vec3(0.0, 0.0, 0.0),
   normal = Vec3(0.0, 1.0, 0.0)
)

#create a wall at the top of the model:
sim.createWall (
   name = "top_wall",
   posn = Vec3(0.0, Ylength, 0.0),
   normal = Vec3(0.00, -1.0, 0.0)
)

#create a wall at the bottom of the model:
sim.createWall (
   name = "left_wall",
   posn = Vec3(0.0, 0.0, 0.0),
   normal = Vec3(1.0, 0.0, 0.0)
)

#create a wall at the bottom of the model:
sim.createWall (
   name = "right_wall",
   posn = Vec3(Xlength, 0.0, 0.0),
   normal = Vec3(-1.0, 0.0, 0.0)
)

#create rotational elastic-brittle bonds between particles:
pp_bonds = sim.createInteractionGroup (
   BrittleBeamPrms(
      name="pp_bonds",
      youngsModulus= YoungMod,
      poissonsRatio= Poisson,
      cohesion= Cohes,
      tanAngle= TanPhi,
      tag=1
   )
)

#initialise frictional interactions for unbonded particles:
sim.createInteractionGroup (
   FrictionPrms(
      name="friction",
      youngsModulus= YoungMod,
      poissonsRatio= Poisson,
      dynamicMu= MuDyn,
      staticMu= MuSta
   )
)

#create an exclusion between bonded and frictional interactions:
sim.createExclusion (
   interactionName1 = "pp_bonds",
   interactionName2 = "friction"
)

#specify elastic repulsion from the bottom wall:
sim.createInteractionGroup (
   NRotBondedWallPrms (
      name = "bottom_wall_repel",
      wallName = "bottom_wall",
      normalK = YoungMod,
      particleTag = 2
   )
)

#specify elastic repulsion from the top wall:
sim.createInteractionGroup (
   NRotBondedWallPrms (
      name = "top_wall_repel",
      wallName = "top_wall",
      normalK = YoungMod,
      particleTag = 3
   )
)

#specify elastic repulsion from the top wall:
sim.createInteractionGroup (
   NRotElasticWallPrms (
      name = "left_wall_repel",
      wallName = "left_wall",
      normalK = YoungMod
   )
)

#specify elastic repulsion from the top wall:
sim.createInteractionGroup (
   NRotElasticWallPrms (
      name = "right_wall_repel",
      wallName = "right_wall",
      normalK = YoungMod
   )
)

#add translational viscous damping:
sim.createInteractionGroup (
   LinDampingPrms(
      name="damping1",
      viscosity=2,
      maxIterations=100
   )
)

#add rotational viscous damping:
sim.createInteractionGroup (
   RotDampingPrms(
      name="damping2",
      viscosity=2,
      maxIterations=100
   )
)

#add a wall loader to move the top wall:
wall_loader1 = WallLoaderRunnable(
   LsmMpi = sim,
   wallName = "right_wall",
   vPlate = Vec3 (-0.1, 0.0, 0.0),
   startTime = TotalTimestep/5,
   rampTime = TotalTimestep/4
)
sim.addPreTimeStepRunnable (wall_loader1)

#add a wall loader to move the bottom wall:
wall_loader2 = WallLoaderRunnable(
   LsmMpi = sim,
   wallName = "left_wall",
   vPlate = Vec3 (0.1, 0.0, 0.0),
   startTime = TotalTimestep/5,
   rampTime = TotalTimestep/4
)
sim.addPreTimeStepRunnable (wall_loader2)

#add ServoWallLoaderRunnables to apply constant normal stress:
servo_loader1 = ServoWallLoaderRunnable(
   LsmMpi = sim,
   interactionName = "top_wall_repel",
   force = Vec3 (0.0,(-0.5*Xlength*ConfiningP), 0.0),
   startTime = 0,
   rampTime = 400000
)
sim.addPreTimeStepRunnable (servo_loader1)

#add ServoWallLoaderRunnables to apply constant normal stress:
servo_loader2 = ServoWallLoaderRunnable(
   LsmMpi = sim,
   interactionName = "bottom_wall_repel",
   force = Vec3 (0.0, (0.5*Xlength*ConfiningP), 0.0),
   startTime = 0,
   rampTime = 400000
)
sim.addPreTimeStepRunnable (servo_loader2)

#create a FieldSaver to store number of bonds:
sim.createFieldSaver (
   InteractionScalarFieldSaverPrms(
      interactionName="pp_bonds",
      fieldName="count",
      fileName="nbonds.dat",
      fileFormat="SUM",
      beginTimeStep=0,
      endTimeStep=TotalTimestep,
      timeStepIncr=1000
   )
)

#create a FieldSaver to store the total kinetic energy of the particles:
sim.createFieldSaver (
   ParticleScalarFieldSaverPrms(
      fieldName="e_kin",
      fileName="ekin.dat",
      fileFormat="SUM",
      beginTimeStep=0,
      endTimeStep=TotalTimestep,
      timeStepIncr=1000
   )
)

#create a FieldSaver to store potential energy stored in bonds:
sim.createFieldSaver (
   InteractionScalarFieldSaverPrms(
      interactionName="pp_bonds",
      fieldName="potential_energy",
      fileName="epot.dat",
      fileFormat="SUM",
      beginTimeStep=0,
      endTimeStep=TotalTimestep,
      timeStepIncr=1000
   )
)

#create a FieldSaver to wall positions:
posn_saver = WallVectorFieldSaverPrms(
   wallName=["bottom_wall", "top_wall", "left_wall", "right_wall"],
   fieldName="Position",
   fileName="out_Position.dat",
   fileFormat="RAW_SERIES",
   beginTimeStep=0,
   endTimeStep=TotalTimestep,
   timeStepIncr=1000
)
sim.createFieldSaver(posn_saver)

#create a FieldSaver to wall forces:
force_saver = WallVectorFieldSaverPrms(
   wallName=["bottom_wall", "top_wall", "left_wall", "right_wall"],
   fieldName="Force",
   fileName="out_Force.dat",
   fileFormat="RAW_SERIES",
   beginTimeStep=0,
   endTimeStep= TotalTimestep,
   timeStepIncr=1000
)
sim.createFieldSaver(force_saver)

#add a CheckPointer to store simulation data:
sim.createCheckPointer (
   CheckPointPrms (
      fileNamePrefix = "snapshot",
      beginTimeStep = 0,
      endTimeStep = TotalTimestep,
      timeStepIncr = 100000
   )
)

#execute the simulation:
sim.run()