Force in a facet

Hello launchpad,

I am new to Yade. I have read the doc, and the examples.
I am trying to simulate a rock cutting and record :
1) the forces acting on the sample, and
2) the cracks

Here is my code:

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division

from yade import plot,pack,timing, utils, geom
import time, sys, os, copy

#import matplotlib
#matplotlib.rc('text',usetex=True)
#matplotlib.rc('text.latex',preamble=r'\usepackage{concrete}\usepackage{euler}')

# default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
 young=24e9, #kn
 poisson=.2, #ks

 sigmaT=3.5e6,
 frictionAngle=atan(0.8),
 epsCrackOnset=1e-4,
 relDuctility=30,

 intRadius=1.5,
 dtSafety=.8,
 damping=0.4,
 #strainRateTension=.05,
 #strainRateCompression=.5,
 #setSpeeds=True,
 # 1=tension, 2=compression (ANDed; 3=both) den xreiazontai, prokeitai gia dokimh kophs
 #doModes=2,

 specimenLength=1,
 sphereRadius=0.99e-3, #0.35e-3

 # isotropic confinement (should be negative)
 isoPrestress=0,
)

from yade.params.table import *

if 'description' in O.tags.keys(): O.tags['id']=O.tags['id']+O.tags['description']
#material properties
concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800,sigmaT=sigmaT,
     relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress))

sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((-75e-3*specimenLength,-27e-3*specimenLength,-6.35e-3*specimenLength),(75e-3*specimenLength,27e-3*specimenLength,6.35e-3*specimenLength)),spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation(material=concreteId)

bx = geom.facetBox(center=(-80e-3,27e-3,0.0),extents=(3e-3,-5e-3,6.35e-3),orientation=Quaternion((0,0,1),-pi/36), wallMask=(2+8))
O.bodies.append(bx)
for facet in bx:
 facet.state.blockedDOFs='xyzXYZ'
 facet.state.vel=(20.0,0,0)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),Bo1_Facet_Aabb()],verletDist=.05*sphereRadius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc'), Ig2_Facet_Sphere_ScGeom()],
  [Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=1)],
  [Law2_ScGeom_CpmPhys_Cpm()],
 ),
 NewtonIntegrator(damping=damping,label='damper'),
 PyRunner(command='addPlotData()',realPeriod=2), #CpmStateUpdater(realPeriod=.5),
]

def addPlotData():
   for facet in bx:
      #Fn = facet.shape.normal
      Fn=O.forces.f(facet.id)[2]
      disp = facet.state.pos[2]
   plot.addData(i=O.iter,Fn=Fn,disp=disp)
plot.plots={'displacement':('Fn',)}
plot.plot()

O.dt=dtSafety*PWaveTimeStep()