stress field

Hi All,

I want to use the TesselationWrapper command to get the micro-stress, it seems that the code works fine. it can export the vtk.file.

however, all the values in each file are the same.

what's more, in my code, I want to save the files every 10000 steps. while the vtk.flies save every step.

here is the MWE:

#############################
from yade import pack, plot
from yade import export
from yade import utils
compFricDegree=45.0
finalFricDegree=19.5
rate=-5
damp=0.7
target_strain = 0.15
stabilityThreshold=0.01
young=4e8
mn,mx=Vector3(0,0,0),Vector3(4e-3,8e-3,4e-3)
O.materials.append(CohFrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),isCohesive=False,alphaKr=0.2,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True,density=2648,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
psdSizes=[0.12e-3,0.14e-3,0.16e-3,0.18e-3,0.19e-3,0.20e-3,0.22e-3,0.24e-3,0.29e-3,0.38e-3]
psdCumm=[0,11.1,22.2,33.3,44.4,55.5,66.6,77.7,88.8,100]
sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,num=5000,distributeMass=1,seed=1)
sp.toSimulation(material='spheres')
O.dt=.5*PWaveTimeStep()
triax=TriaxialStressController(
 maxMultiplier=1.0+2e5/young,
 finalMaxMultiplier=1.0+2e4/young,
 thickness=0,
 stressMask=7,
 internalCompaction=False
)
newton=NewtonIntegrator(damping=damp)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
   [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
   [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True,label='cohesiveLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8),
 triax,
 newton,
]
triax.goal1=triax.goal2=triax.goal3=-100000
while 1:
 O.run(3000,True)
 unb=unbalancedForce()
 print 'unbalanced force: ',unb,'mean stress: ',triax.meanStress
 if unb<stabilityThreshold and abs(-100000-triax.meanStress)/-100000<0.01:
  break
triax.wall_bottom_activated=True
triax.wall_top_activated=True
triax.wall_left_activated=True
triax.wall_right_activated=True
triax.wall_back_activated=True
triax.wall_front_activated=True
setContactFriction(radians(finalFricDegree))
re11=-triax.strain[0]
re22=-triax.strain[1]
re33=-triax.strain[2]
triax.stressMask = 5
newton=NewtonIntegrator(damping=0.0)
triax.goal2=rate
triax.strainRate2=5
triax.strainRate1=triax.strainRate3=1000.0
def stop_loading():
 if -triax.strain[1]-re22> target_strain:
  O.pause()

def stop_loading():
 if -triax.strain[1]-re22> target_strain:
  O.pause()
################## here is the stess field part ###############
TW=TesselationWrapper()
TW.setState()
TW.computeVolumes()
s=bodyStressTensors()
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.mystress = s[b.id]
vtk = export.VTKExporter("stress-needed")

step1 = 0
while -triax.strain[1]-re22< target_strain:
 step1 = step1+10000 ##### I want to save it every 10000 steps.
 vtk.exportSpheres(ids='all',what=[('stress','b.mystress')],useRef=True)

O.engines =O.engines[0:6]+[PyRunner(command='stop_loading()',iterPeriod=10000)]+O.engines[6:7]

O.run(500000,True)

##########################

thanks,
yong