Export Deviatoric loading phase data

Hi, everybody!
I am a Yade newbie and now I want to export the data from each step of the deviatoric loading phase. Here is my code. What should I do? Thank you for your help.

from yade import pack,plot,qt
import math
############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
        num_spheres=1000,# number of spheres
        compFricDegree = 30, # contact friction during the confining phase
        key='_triax_base_', # put you simulation's name here
        unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key
targetPorosity = 0.43 #the porosity we want for the packing
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=-0.02 # loading rate (strain rate)
damp=0.2 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=5e6 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing
## create materials for spheres and plates
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
wallIds=O.bodies.append(walls)
psdSizes,psdCumm =[0.05,0.1],[0.5,1]# [0.02,0.04,0.045,0.05,0.06,0.08,0.12],[0.,0.1,0.3,0.3,0.3,0.7,1.]
sp = pack.SpherePack()
sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True)
#Add material
sp.toSimulation(material='spheres')
triax=TriaxialStressController(
        ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
        ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
        maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
        finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
        thickness = 0,
        ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
        ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
        ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
        ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
        ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
        stressMask = 7,
        internalCompaction=False, # If true the confining pressure is generated by growing particles
)
newton=NewtonIntegrator(damping=damp)
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        triax,
        TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
        newton,
        PyRunner(command='addPlotData()',realPeriod=2),
  # VTKRecorder(fileName='/home/ma/vtk/gravity-',recorders=['all'],label='vvtk',iterPeriod=100)
]
def addPlotData():
  plot.addData(i=O.iter,
                     e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
                    s11=-triax.stress(triax.wall_right_id)[0],
                    s22=-triax.stress(triax.wall_top_id)[1],
                    s33=-triax.stress(triax.wall_front_id)[2],
                    Ek=utils.kineticEnergy(),coordNum=utils.avgNumInteractions(),unForce=utils.unbalancedForce())
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()
triax.goal1=triax.goal2=triax.goal3=-10000
while 1:
  O.run(1000, True)
  ##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
    break
print "### Isotropic state saved ###"
#O.run(1000,True)
#O.save('test2.yade')
triax.internalCompaction=False
## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))
##set stress control on y and z, we will impose strain rate on x
triax.stressMask = 5
##now goal2 is the target strain rate

triax.goal1=-10000
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal2=rate
triax.goal3=-10000
newton.damping=0.1
#O.run(100,True)
O.saveTmp()
##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
plot.plots={'e22':('s11','s22','s33',None,'ev')}
plot.plot()
O.run(1000,True)
for b in O.bodies:
      print "Body",b.id
      print " displ",b.state.displ()
      print " vel ",b.state.vel
      for i in b.intrs():
        print " Intr with",i.id2 if i.id1 == b.id else i.id1
        print " cp",i.geom.contactPoint
        print " f ",i.phys.normalForce+i.phys.shearForce