Particle breakage of clump using JCFpm

Hi! I want to model the particle breakage of clumps during the one-dimensional compression in a cylindrical boundary. My steps are as followed:
Step 1: I import 800 GTS files into Yade and generated one-on-one 800 clumps. That is to say that each GTS corresponds to a unique clump.
Step 2: I create a cylindrical boundary to cover those clumps and provide two plates (bottom one fixed, top one applied with force). Note that the force applied on the top plate is a multi-load, which varies with the strain of the sample.
Step 3: I assign JCFpm material for each clump and material of steel for the boundary and plates.
Step 4: record what we want to know and implement simulation.

Here, my questions are as follows:
1) Is my whole procedure appropriate to simulate the one-dimensional compression while considering particle breakage?
2) For Step 2, I would like to know how to assign multi-loads is correct. It seems currently the top plate is not quasi-static during loading. Do I need to clear the speed for every number of calculation steps?
3) For the whole simulation, how do we know how many clumps are broken and can we record the isolated clumps at the final?

The MWM is presented below (I don't know how to provide the gts files without an external link, so maybe randomly generating some clumps can replace them.):

from __future__ import print_function
from yade import pack, ymport
import gts, os.path, locale, plot, random
import numpy as np

locale.setlocale(
        locale.LC_ALL, 'en_US.UTF-8'
) #gts is locale-dependend. If, for example, german locale is used, gts.read()-function does not import floats normally

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

# The following lines are used parameter definitions
readParamsFromTable(
 # directory
 dir_gts = '/home/kyle/Yade/install/bin/gts', # directory of storing gts files
 dir_vtk = '/home/kyle/Yade/install/bin/vtk', # directory of storing vtk files
 dir_txt = '/home/kyle/Yade/install/bin/txt', # directory of storing txt files
 fileType = '.gts', # type of file

        # type of generating spheres ['hexa','dense']
        genType = 'hexa',

        # material parameters
        young_w = 210e9, # Young's modulus of plates and walls
        young_g = 180e6, # Young's modulus of grains [1]
        den_w = 7874, # density of plates and walls
        den_g = 1162, # density of grains
        poi_w = 0.3, # possion's ratio of plates and walls
        poi_g = 0.25, # possion's ratio of grains
        friAngle_w = 15, # friction angle of plates and walls
        friAngle_g = 19, # friction angle of grains

        # Parameters of cylindrical walls
        x_cyl = 0.0106, # x-coordinate of center of cylindrical walls
        y_cyl = 0.01005, # y-coordinate of center of cylindrical walls
        z_cyl = 0, # z-coordinate of center of cylindrical walls
        r_cyl = 0.0097, # radius of the cylinder
        h_cyl = 0.04151, # height of the cylinder
        n_cyl = 100, # number of boxes forming cylinder
        t_cyl = 0.0001, # thickness of the boxes

        #Parameters of lower and upper plates (square box)
        w_plate = 0.015, # width of plates
        t_plate = 0.0001, # thickness of plates

        # applied stress
        stress = [12500, 25000, 50000, 100000, 200000, 300000, 400000, 800000, 1600000, 3200000], # target applied stress

        # target displacement
        displacement = [3.6e-5, 6.8e-5, 0.00014, 0.00028, 0.00047, 0.000642, 0.000796, 0.001342, 0.002441, 0.003989] # target displacement

)
from yade.params.table import *

# create materials for spheres and walls
wallMat = O.materials.append(FrictMat(young = young_w, poisson = poi_w, frictionAngle = friAngle_w, density = den_w, label = 'walls'))
#sphereMat = O.materials.append(FrictMat(young = young_g, poisson = poi_g, frictionAngle = friAngle_g, density = den_g, label = 'spheres'))

def sphereMat():
 return JCFpmMat(
         type=1,
         young = young_g,
         frictionAngle = radians(friAngle_g),
         density = den_g,
         poisson = poi_g,
         tensileStrength = 1e6,
         cohesion=1e6,
         jointNormalStiffness=1e7,
         jointShearStiffness=1e7,
         jointCohesion=1e6,
         jointFrictionAngle=radians(20),
         jointDilationAngle=radians(0)
 ) ## Rq: density needs to be adapted as porosity of real rock is different to granular assembly due to difference in porosity (utils.sumForces(baseBodies,(0,1,0))/(Z*X) should be equal to Gamma*g*h with h=Y, g=9.82 and Gamma=2700 kg/m3

###############################################
### IMPORTING GRAINS AND CREATING CLUMPS ###
###############################################

##### a function that searching all .gts files
def SearchFiles(directory, ftype):
 fileList = []
 for root, subDirs, files in os.walk(directory):
  for fileName in files:
   if fileName.endswith(ftype):
    fileList.append(os.path.join(root, fileName))
 return fileList
#####

# import gts and obatin file list
fileList = SearchFiles(dir_gts,fileType)

##### a function that determining the radius of spheres
def createSphere(pred, ndivmin, spacing, genType):
 #dim = pred.dim()
 #minDim = min(dim[0], dim[1], dim[2])
 aabb = pred.aabb()
 dim0 = aabb[1][0] - aabb[0][0]
 bodyList = []
 ndiv = ndivmin
 if genType == 'hexa':
  while len(bodyList)==0:
   radius = dim0 / ndiv # get some characteristic dimension, use it for radius
   bodyList = O.bodies.append(pack.regularHexa(pred, radius = radius, gap = -radius/4, color = [random.random(),random.random(),random.random()], material = sphereMat))
   ndiv = ndiv + spacing
 else:
  while len(bodyList)==0:
   radius = dim0 / ndiv # get some characteristic dimension, use it for radius
   sp = pack.randomDensePack(pred,radius=radius,returnSpherePack=True, color = [random.random(),random.random(),random.random()], material = sphereMat)
   bodyList = sp.toSimulation()
   ndiv = ndiv + spacing
 return bodyList
#####

# import gts and create clump for each .gts file
sphereList = []
for gtsName in fileList:
 surf = gts.read(open(gtsName))
 if surf.is_closed():
  print(gtsName)
  pred = pack.inGtsSurface(surf)
  bodyList = createSphere(pred, 20, 2, genType)
  idClump = O.bodies.clump(bodyList)
 del surf
 del pred
 del bodyList
O.bodies.updateClumpProperties()
print('Number of bodies:', len(O.bodies))

###################################################################
### CREATING CYLINDRICAL WALLS, BOTTOM PLATE, AND TOP PLATE ###
###################################################################

##### a function that creating cylindrical walls
def cylSurf(center,radius,height,nSegments=12,thick=0,**kw):
   """creates cylinder made of boxes. Axis is parallel with z+
   center ... center of bottom disc
   radius ... radius of cylinder
   height ... height of cylinder
   nSegments ... number of boxes forming cylinder
   thick ... thickness of the boxes
   **kw ... keywords passed to box function (e.g. material, color ...)
   """
   # just converts (a,b,c) to Vector3(a,b,c) to enable "center + ..."
   center = Vector3(center)
   # nSegments angles to define individual boxes, from 0 to 2*pi
   angles = [i*2*pi/nSegments for i in range(nSegments)]
   # centers of boxes. Vector 3(radius,0,0) rotated by angle 'a'
   # z coordinate is .5*height
   # center + shift = center of corresponding box
   pts = [center + Vector3(radius*cos(a),radius*sin(a),.5*height) for a in angles] # centers of plates
   ret = []
   for i,pt in enumerate(pts): # for each box center create corresponding box
      # "length" of the box along circumference, cylinder circumference / nSegments
      l = pi*radius/nSegments
      # extents, half dimensions along x,y,z axis
      es = (.5*thick,l,.5*height)
      # by default box is axis-aligned, we need to rotate it along z axis (0,0,1) by angle i*2*pi/nSegments
      ori = Quaternion((0,0,1),i*2*pi/nSegments)
      # have a look at utils.box documentation. pt=center of box, es=extents (half dimension), ori=orientation
      ret.append(box(pt,es,ori,**kw))
   return ret
#####

# create cylindrical walls and fix them
surf = cylSurf([x_cyl, y_cyl, z_cyl], r_cyl, h_cyl, nSegments = n_cyl, thick = t_cyl, material = wallMat, wire=True)
O.bodies.append(surf)
for b in surf: b.state.blockedDOFs = 'xyzXYZ'
for b in surf: b.state.vel = (0, 0, 0)

# create the bottom plate and fix it
bottom_plate = O.bodies.append(box([x_cyl, y_cyl, z_cyl], (w_plate, w_plate, t_plate), material = wallMat))
O.bodies[bottom_plate].state.vel = (0, 0, 0)
O.bodies[bottom_plate].state.blockedDOFs = 'xyzXYZ'

# create the top plate
top_plate=O.bodies.append(box([x_cyl, y_cyl, h_cyl], (w_plate, w_plate, t_plate), material = wallMat))

# apply initial force
force = np.dot(stress, w_plate * w_plate) # corresponding applied force
global deltaF
deltaF = force[0]
O.forces.setPermF(top_plate,(0, 0, -deltaF))

# target displacement
global numIter
numIter = 0
global maxIter
maxIter = 9
global dispT
dispT = displacement[numIter]
global disp
disp = 0
global dispE
dispE = displacement[-1]

# calculate initial void ratio
vol0 = h_cyl * pi * r_cyl * r_cyl
volG0 = sum([b.state.mass for b in O.bodies if b.isClump])/den_g
volV0 = vol0 - volG0
e0 = volV0/volG0
print(volG0, e0)
print('numIter', 'dispT', 'disp', 'disp/dispT', 'velTopPlate')

############################
### DEFINING ENGINES ###
############################

O.engines = [
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1),
    Bo1_Box_Aabb()]),
 InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1),
                 Ig2_Box_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys(),
                 Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
                [Law2_ScGeom_FrictPhys_CundallStrack(),
                 Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True)]),
         GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.8),
 PyRunner(iterPeriod = 10,command = "checkDisplacement(displacement, force, h_cyl)"),
 VTKRecorder(iterPeriod = 500, recorders = ['jcfpm', 'cracks', 'spheres', 'boxes', 'stress', 'id', 'clumpId', 'coordNumber', 'colors'], fileName = dir_vtk + '/'),
 NewtonIntegrator(gravity = [0, 0, -9.81], damping = 0.4)
]
O.dt= 0.1 * PWaveTimeStep()

# check whether the target displacement has already met. if so, update applied stress
def checkDisplacement(displacement, force, h_cyl):
 global numIter
 global maxIter
 global disp
 global dispT
 global dispE
 global deltaF
 # at the very start, applied stress is the first of value of array 'stress'
 t = O.time
 H = O.bodies[top_plate].state.pos[2]
 disp = h_cyl - H
 velTopPlate = O.bodies[top_plate].state.vel[2]
 if disp > dispT:
  if numIter < maxIter:
   numIter = numIter + 1
   dispT = displacement[numIter]
   #deltaF = force[numIter] - force[numIter - 1]
   O.forces.setPermF(top_plate,(0, 0, -force[numIter]))
  else:
   print(numIter, dispT, disp, disp/dispT, velTopPlate, O.forces.permF(top_plate)[2], O.forces.f(top_plate)[2])
   f = open(dir_txt + '/data.txt', 'a')
   f.write('{} {} {} {} {} {} {}\n'.format(O.iter, numIter, dispT, disp, velTopPlate, O.forces.permF(top_plate)[2], H))
   f.close()
   plot.saveDataTxt(O.tags['d.id'] + '.txt')
   O.pause()
 print(numIter, dispT, disp, disp/dispT, velTopPlate, O.forces.permF(top_plate)[2], O.forces.f(top_plate)[2])
 plot.addData(t=t,disp=disp)
 f = open(dir_txt + '/data.txt', 'a')
 f.write('{} {} {} {} {} {} {}\n'.format(O.iter, numIter, dispT, disp, velTopPlate, O.forces.permF(top_plate)[2], H))
 f.close()

##########################
### RUN SIMULATION ###
##########################
plot.plots = {'t':('disp')}
plot.plot()
O.run()