Is possible to plot the unbreakable bonds during UTS in paraview

Hello to everybody, i hope everyone is fine during Convid!!

To be honest, i beleive the title of my question says everything.
I am conducting UTS simulations and i want to process the reuslts in paraview.
I dont know how to plot the unbreakable zone of bonds in paraview?
Is any way fir something like that to be achieved?

Thank you.
I am showing the code below...

# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-
from yade import ymport, plot

#### Simulation of a uniaxial test (compression or tension depending on the sign of the loading rate "rate" defined below)
# everything is documented at https://yade-dem.org/doc/index-toctree.html -> use the quick serach bar on the left of the page

################# SIMULATIONS DEFINED HERE

#### Simulation Control
rate=0.01 # deformation rate (0.003 for tension, -0.03 for compression)
iterMax=10000 # maximum number of iterations
saveVTK=4000 # saving output files for paraview
OUT='tensionTest_r0.01' # output files

intR=1.00 # allows near neighbour interaction (can be adjusted for every packing)
DENS=5265 # 2640 toy petrwmatos# could be adapted to match material density: dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> packing porosity as to be computed?
YOUNG=6e2
FRICT=13
ALPHA=2.5
TENS=2e6
COH=1.5*TENS

#### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
D=5e-2
H = 2.5*D

pred=pack.inCylinder((0,0,0),(0,H,0),D)

sps=SpherePack()
sp=pack.randomDensePack(pred,spheresInCell=10000,radius=1e-3,rRelFuzz=.34,memoizeDb='/tmp/test.sqlite',returnSpherePack=True)
sp.toSimulation(material=sphereMat) # up to you to define another sample here, e.g., with randomDensePack or anything else.

#### Particle Characteristics
R=0
Rmax=0
nbSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   nbSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
Rmean=R/nbSpheres
print "\n"
print ('nbSpheres=',nbSpheres,' | Rmean=',Rmean,'| Rmax=',Rmax,)
print "\n"

#### Set up boundary conditions (see https://yade-dem.org/doc/yade.utils.html?highlight=uniaxialtestfeatures#yade.utils.uniaxialTestFeatures)
bb=utils.uniaxialTestFeatures()
negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

################# ENGINES DEFINED HERE

O.engines=[
 ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom')],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
 ),
 UniaxialStrainer(strainRate=rate,axis=longerAxis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=0,stopStrain=0.1,dead=1,label='strainer'),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5, defaultDt=utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.4,label='newton'),
 PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
 PyRunner(iterPeriod=int(100), command='Stop()',label='data'),
        VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
]

################# RECORDER DEFINED HERE

def recorder():
    yade.plot.addData({'i':O.iter,
         'eps':strainer.strain,
         'sigma':strainer.avgStress,
         'tc':interactionLaw.nbTensCracks,
         'sc':interactionLaw.nbShearCracks,
         'te':interactionLaw.totalTensCracksE,
         'se':interactionLaw.totalShearCracksE,
         'unbF':utils.unbalancedForce()})
    plot.saveDataTxt(OUT)

def Stop():
    sigma=plot.data['sigma']
    if abs(strainer.avgStress) <= 0.75*max(sigma):
        O.pause()

# if you want to plot during simulation
plot.plots={'eps':('sigma')}
plot.plot()

################# PREPROCESSING
print " Number of particles: %.2f" %(len(O.bodies))
print " Porosity of sample is: %.4f" %(utils.porosity())
print "\n"

#### manage interaction detection factor during the first timestep and then set default interaction range ((cf. A DEM model for soft and hard rock, Scholtes & Donze, JMPS 2013))
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

# reinforce bonds at boundaries for tensile tests
dim=utils.aabbExtrema()
if rate>0:
 layerSize=0.15
 for o in O.bodies:
  if isinstance(o.shape,Sphere):
   if ( o.state.pos[longerAxis]<(dim[0][longerAxis]+layerSize*(dim[1][longerAxis]-dim[0][longerAxis])) ) or ( o.state.pos[longerAxis]>(dim[1][longerAxis]-layerSize*(dim[1][longerAxis]-dim[0][longerAxis])) ) :
    o.shape.color=(1,1,1)
 for i in O.interactions:
  if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
   if O.bodies[i.id1].shape.color==(1,1,1) or O.bodies[i.id2].shape.color==(1,1,1) :
    i.phys.FnMax*=100
    i.phys.FsMax*=100

#### coordination number after the creation of contacts at first step. (Before = 0)
print " Number of Interactions per Particle: %.2f" %(utils.avgNumInteractions())#CoordinationNumber
print "\n"
#### coordination number verification
numSSlinks=0
numCohesivelinks=0
for i in O.interactions:
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
      numSSlinks+=1
    if i.phys.isCohesive :
      numCohesivelinks+=1

print ("nblinks=", numSSlinks, " | nbCohesivelinks=", numCohesivelinks, " || Kcohesive=", 2.0*numCohesivelinks/nbSpheres)
print "\n"
################# SIMULATION REALLY STARTS HERE
strainer.dead=0
#O.run(iterMax)