flow.nCells() != "cell number" recorded by flow.saveVtk()

Hello,

One question for FlowEngine and DFNFlowEngine: why the flow.nCells() != "cell number" recorded by flow.saveVtk()?

I would like to record the Ids of initial fractured cells and then impose different pressures within them, but I find that the ids (and vertices) are not consistent given by these two ways.

Below is an example.
Thank you in advance for your help.
Best regards,
Lingran

from yade import pack, ymport

### material
def sphereMat(): return JCFpmMat(type=1,density=3000,young=9e9,poisson=0.2,frictionAngle=radians(18))
def wallMat(): return JCFpmMat(type=0,density=3000,young=9e9,poisson=0.2,frictionAngle=radians(0))

### bodies
mn,mx=Vector3(0,0,0),Vector3(1,1,1)
O.bodies.append(aabbWalls([mn,mx],oversizeFactor=1.5,thickness=0.1,material=wallMat))
sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((0.,0.,0.),(1,1.,1.)),spheresInCell=100,radius=1/10.,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation(material=sphereMat)

dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

### engines
flow=FlowEngine(
        isActivated=1,
        ### choose solver to use (0: Gauss Seidel, 1: Taucs, 2: Pardiso, 3: CHOLMOD)
        useSolver=3, # 3 should be used by default but does not work everytime...
         boundaryUseMaxMin = [0,0,0,0,0,0],
        bndCondIsPressure = [1,1,0,0,0,0],
        bndCondValue = [1,0,0,0,0,0],
        permeabilityFactor=1,
        viscosity=0.001,
        #debug=1
)

intR=1.1
O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')])
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
 )
 ,flow
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8)
 ,NewtonIntegrator()
]

### simulation starts here

## to block the particles (better for permeability test)
for o in O.bodies:
 o.state.blockedDOFs+='xyzXYZ'

O.run(1,True)
## getBoundaryFlux get the total discharge [m3/s]
Qin = flow.getBoundaryFlux(0)
Qout = flow.getBoundaryFlux(1)

permeability = abs(Qout)*flow.viscosity*X/(Y*Z) # !!! if Pout=1, Pin=0
permeability2 = flow.averageVelocity()*flow.viscosity*X # !!! if Pout=1, Pin=0
conductivity = permeability*1000*9.82/flow.viscosity # K=rho*g*k/nu
print "Qin=",Qin," Qout=",Qout," ARE THEY EQUAL? IF NOT-> NO FLOW!"
print "Permeability [m2]=",permeability," || Hydraulic conductivity [m/s]=",conductivity
flow.saveVtk() # to see the result in Paraview

for i in range(flow.nCells()):
   print i, flow.getVertices(i)