Some fundamental questions about the Brazil Splitting Test example file.

Hello everyone, I have some fundamental questions regarding the Brazil Splitting Test example file[1]. I made some minor modifications on the basis of the source code (without altering the main code). Although the code runs, I have some uncertainties regarding the details. I would greatly appreciate any assistance in addressing these questions.

1. Regarding the cohesive state of particles in the CPM material, is it accurate to say that all contacting particles form cohesive connections? In the engine section, the aabbEnlargeFactor and interactionDetectionFactor are set to 1.5. Does this imply that when the distance between the centers of two particles with CPM material is less than 1.5 times the sum of their radii, a cohesive connection, also referred to as a bond, is established? If this understanding of cohesion is correct, what does the code annotated as "now reset the interaction radius and go ahead" mean in the subsequent part of the code? Why is the factor changed to 1? I'm seeking clarification on the meaning of this comment and code.
2. When configuring the material, I specified that once a cohesive connection is broken, it cannot be reestablished. Therefore, I would like to know the time and location of each cohesive bond breakage. How can I obtain this information? Additionally, how can I visualize the positions of these cracks using VTK? I attempted to use parameters like currentinteraction or currentCohesiveBond to determine the number of cracks, but it seems ineffective. Firstly, I assumed their values should be consistent, but in reality, they are not entirely the same. Secondly, as the simulation progresses, they do not seem to decrease and sometimes even increase, which appears illogical.
3. In this model, all initial particles are in blue, and as loading progresses, some particles turn purple-red. What do these purple-red particles represent? Do they indicate the cumulative energy loss of all connecting bonds around them? It's worth noting that in some other simulations, spherical particles are colored and do not show changes in force. Is there a command controlling this variation?
4. When the simulation load reaches a certain point, the stopIfDamaged function halts the loading. However, when I remove this function, the model continues to load, and the specimen quickly explodes, deviating from real experimental conditions. I believe the specimen should fracture into several pieces, each composed of many particles, rather than dispersing almost evenly into particles with minimal connections. Additionally, in the Brazil Splitting Test, the fracture should initially appear at the center of the specimen. In the simulation, however, most initial fractures occur at the loading end. How should I modify the code to make the simulation more realistic?

Thank you in advance!

Below is the code:
-------------------------------------------------------------------------
# -*- encoding=utf-8 -*-
from __future__ import print_function
from yade import plot, pack
import xlsxwriter
"""
A simple script of a Brazilian splitting test.
A cylinder with x axis form 0 to specimenLength.
Load in z direction by z-perpendicular walls.

Code strongly inspired by uniax.py
"""

# default parameters or from table
readParamsFromTable(
        noTableOk=True, # unknownOk=True,
        young=24e9,
        poisson=.2,
        sigmaT=3.5e6,
        frictionAngle=atan(0.8),
        epsCrackOnset=1e-4,
        relDuctility=30,
        intRadius=1.5,
        dtSafety=.8,
        strainRate=1,
        specimenLength=.05,
        specimenRadius=.05,
        sphereRadius=2.5e-3,
)
from yade.params.table import *

# material
concreteId = O.materials.append(
        CpmMat(
                young=young,
                poisson=poisson,
                frictionAngle=frictionAngle,
                epsCrackOnset=epsCrackOnset,
                relDuctility=relDuctility,
                sigmaT=sigmaT,
        )
)

# spheres
sp = pack.randomDensePack(
        pack.inCylinder((0, 0, 0), (specimenLength, 0, 0), specimenRadius),
        radius=sphereRadius,
        spheresInCell=500,
        memoizeDb='/tmp/packing-brazilian.db',
        returnSpherePack=True
)
sp.toSimulation()

# walls
zMin, zMax = [pt[2] for pt in aabbExtrema()]
wallIDs = O.bodies.append([wall((0, 0, z), 2) for z in (zMin, zMax)])
walls = wallMin, wallMax = [O.bodies[i] for i in wallIDs]
v = strainRate * 2 * specimenRadius
wallMin.state.vel = (0, 0, +v)
wallMax.state.vel = (0, 0, -v)

timeStr = time.strftime("%Y-%m-%d %X", time.localtime())
outputDir = 'out_' + '_' + timeStr
nestedDir = 'VTK_'
if not os.path.exists(outputDir):
 os.mkdir(outputDir)
nestedDirPath = os.path.join(outputDir,nestedDir)
if not os.path.exists(nestedDirPath):
 os.mkdir(nestedDirPath)

# engines
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius, label='is2aabb'),
                               Bo1_Wall_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius, label='ss2sc'),
                 Ig2_Wall_Sphere_ScGeom()],
                [Ip2_CpmMat_CpmMat_CpmPhys()],
                [Law2_ScGeom_CpmPhys_Cpm()],
        ),
        NewtonIntegrator(damping=.2),
        CpmStateUpdater(realPeriod=.5),
        VTKRecorder(fileName= outputDir + '/' + nestedDir + '/Brazil-',recorders=['all'],iterPeriod=100),

        PyRunner(iterPeriod=100, command='addPlotData()', initRun=True),
        PyRunner(iterPeriod=100, command='stopIfDamaged()'),
        PyRunner(command='cohesivenum()',iterPeriod=100),
        PyRunner(command='writeData()',iterPeriod=100),
]

# stop condition
def stopIfDamaged():
 if O.iter < 1000: # do nothing at the beginning
  return
 fMax = max(plot.data["f"])
 f = plot.data["f"][-1]
 if f / fMax < .6:
  print("Damaged, stopping.")
  print("ft = ", max(plot.data["stress"]))
  O.pause()

# plot stuff
def addPlotData():
 global t,i,dspl,stress
 # forces of walls. f1 is "down", f2 is "up" (f1 needs to be negated for evlauation)
 f1, f2 = [O.forces.f(i)[2] for i in wallIDs]
 f1 *= -1
 # average force
 f = .5 * (f1 + f2)
 # displacement (2 times each wall)
 wall = O.bodies[wallIDs[0]]
 dspl = 2 * wall.state.displ()[2]
 # stress (according to standard brazilian test evaluation formula)
 stress = f / (pi * specimenRadius * specimenLength)
 # store values
 yade.plot.addData(
         t=O.time,
         i=O.iter,
         dspl=dspl,
         f1=f1,
         f2=f2,
         f=f,
         stress=stress,
 )

##########################
##########################
ballnum=0
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  ballnum=ballnum+1
print('ballnum=',ballnum)

#currentCohesiveBond=0
#currentinteraction=0

def cohesivenum():
 global currentCohesiveBond,currentinteraction
 currentCohesiveBond=0
 currentinteraction=0
 for i in O.interactions:
  if hasattr(i.phys,'isCohesive'):
   currentCohesiveBond+=1
  if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
   currentinteraction+=1
 print('currentinteraction',currentinteraction,end=' ')
 print('currentCohesiveBond',currentCohesiveBond,end=' ')
 print('iter=',O.iter,end=' ')

# save data into *xlsx
t_list = []
dspl_list = []
stress_list = []
currentCohesiveBond_list = []
currentinteraction_list = []
def writeData():
 # write data into lists
 t_list.append(O.time)
 dspl_list.append(dspl)
 stress_list.append(stress)
 currentCohesiveBond_list.append(currentCohesiveBond)
 currentinteraction_list.append(currentinteraction)
 # write data into *xlsx
 workbook = xlsxwriter.Workbook(os.path.join(outputDir, 'Output'+'.xlsx'))
 worksheet = workbook.add_worksheet()
 format = workbook.add_format()
 format.set_align('center')
 worksheet.set_column('A:E',20)
 headings = ['t','dspl','stress','currentCohesiveBond', 'currentinteraction']
 worksheet.write_row('A1',headings,format)
 worksheet.write_column('A2',t_list,format)
 worksheet.write_column('B2',dspl_list,format)
 worksheet.write_column('C2',stress_list,format)
 worksheet.write_column('D2',currentCohesiveBond_list,format)
 worksheet.write_column('E2',currentinteraction_list,format)
 workbook.close()
##########################
##########################

# plot dspl on x axis, stress on y1 axis and f,f1,f2 in y2 axis
plot.plots = {'dspl': ('stress', None, 'f1', 'f2', 'f')}

O.dt = 0.
O.step()
# to create initial contacts
# now reset the interaction radius and go ahead
ss2sc.interactionDetectionFactor = 1.
is2aabb.aabbEnlargeFactor = 1.

# time step
O.dt = dtSafety * PWaveTimeStep()

# run simulation
plot.plot()
#O.run()
-------------------------------------------------------------------------

Thank you all,
Antonia

[1]https://gitlab.com/yade-dev/trunk/-/blob/master/examples/concrete/brazilian.py?ref_type=heads