problem of the command "O.cell.velGrad"

Hi

Firstly, A RVE (goal stress-strain stress(0,0)=0,strain(1,1)=-0.008,stress(0,1)=0, strain(2,:)=0) is calibrated using the following commands:
-----------------------------------------------

#!/usr/bin/python # This is server.py file
import string
from yade import plot,qt
from yade.pack import *
from yade import pack

O.materials.append(CpmMat(damLaw=0,young=150e9,density=4800,frictionAngle=atan(0.8),poisson=0.4,sigmaT=20e8,epsCrackOnset=0.8e-4,relDuctility=300))
initSize=0.028
sp=pack.randomPeriPack(radius=.000625,initSize=Vector3(initSize,initSize,initSize))
O.periodic=True
sp.toSimulation()

O.save('model.gz')

O.dt=1e-7

def plotAddData():
 with open('stressStrain.txt','a') as f:
  f.write(str(p3d.strain[0])+' '+str(p3d.strain[1])+ ' ' + str(p3d.strain[5])+' '+str(p3d.stress[0])+' '+str(p3d.stress[1])+ ' ' + str(p3d.stress[5]))
  f.write('\n')

EnlargeFactor=1.5
EnlargeFactor=1.0
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=EnlargeFactor,label='bo1s')]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=EnlargeFactor,label='ig2ss')],
  [Ip2_CpmMat_CpmMat_CpmPhys()],[Law2_ScGeom_CpmPhys_Cpm()]),
 NewtonIntegrator(),
 Peri3dController(

       nSteps=10000, # how many time steps the simulation will last
       # after reaching nSteps do doneHook action
       doneHook='print "Simulation with Peri3dController finished."; O.pause()',
       label='p3d'
       ),
 PyRunner(command='plotAddData()',iterPeriod=1),
]

p3d.stressMask=0b100001 # prescribed sx,ey,sz,syz,ezx,sxy; e..strain; s..stress
p3d.goal=(0,-0.008,0,0,0,0) #xx, yy, zz, yz, zx, xy
O.step()
bo1s.aabbEnlargeFactor=ig2ss.interactionDetectionFactor=1.
O.run(); O.wait()
------------------------------------------------------

As you can see, the strain and stress at each step are written in the file "stressStrain.txt".

Now I want to use the strain in the "stressStrain.txt" to calculate the stress again using the following commands:

-----------------------------------------------------

#!/usr/bin/python # This is server.py file
import string

from yade.pack import *
from yade import pack
import numpy as np # initial the math of matrix compute (needed in ns)

# initial the global , variables
arrcc =[0]*3 # last step strain(commit)
arrtt =[0]*3 # last step strain(trail)
scc=[0]*3 # last step stress
orst=1 # initialize the orginal save document(trial)
check=[0]*1 # I need to calculate the first step when ag=0
check[0]=0
pertstrain=[0]*3
count=[0]

# load the identical package
O.load('./model.gz') # only load the package, because the O.cell.velGrad is not compatible with the peri3dcontroller eigine
#O.periodic=True ### note!!! add
O.dt=1e-8

ns=100
# see the example of peri3dcontroller, 1.5 is used to speed up the computation
EnlargeFactor=1.5
EnlargeFactor=1.0
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=EnlargeFactor,label='bo1s')],allowBiggerThanPeriod=True),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=EnlargeFactor,label='ig2ss')],
  [Ip2_CpmMat_CpmMat_CpmPhys()],[Law2_ScGeom_CpmPhys_Cpm()]),
 NewtonIntegrator(),

]

O.cell.velGrad=utils.Matrix3(0,0,0,0,0,0,0,0,0)
#O.step()
O.run(1,True)
bo1s.aabbEnlargeFactor=ig2ss.interactionDetectionFactor=1. # O.step and change 1.5 back to 1. , because the concrete model is developed in the 1. environment

def response(msg):
 # change the form of message
 aa = msg.split(',')
 global arrcc
 global arrtt
 global scc
 global orst

     arr=[0]*3 # change form, the current strain

 if orst == 2:
  O.loadTmp('temp')
 orst == 2

     for i in range(0,3):
         arr[i] = string.atof(aa[i])
 ag = [0]*3 # current - last ie. the changed strain
 ag[0] = arr[0] - arrcc[0]
 ag[1] = arr[1] - arrcc[1]
  ag[2] = (arr[2] - arrcc[2]) # incremental strain

 arrcc[0]=arr[0] # update strain
 arrcc[1]=arr[1]
 arrcc[2]=arr[2]

 dstrain = utils.Matrix3(ag[0],ag[2],0,ag[2],ag[1],0,0,0,0) # the goal strain xx xy xz, yx yy yz, zx zy zz
 O.cell.velGrad=dstrain/(ns*O.dt)
 O.run(ns,True)

 stressValue=utils.getStress()
 stressxx=stressValue[0,0] # the return stress xx
 stressyy=stressValue[1,1] # the return stress yy
 stressxy=(stressValue[0,1]+stressValue[1,0])/2 # the return stress xy = (xy + yx)/2
 print stressxx,stressyy,stressxy

 O.saveTmp('temp') # different gauss point has different number(first first ...)

with open('stressStrain.txt','r') as f:
 lines=f.readlines()
flen=len(lines)
for i in lines:
 numbers=i.split()
 strain=map(float,numbers)
 msg = str(strain[0])+","+str(strain[1])+","+str(strain[2])
 response(msg)

--------------------------------------------------------------------------------------------------------

I thought the stress calculated are the same ones in the "stressStrain.txt" (or acceptable small errors ), however, they are not, and the results are seriously affected by the time step O.dt.
I think using the command "O.cell.velGrad=dstrain/(ns*O.dt) "can reach the strain I want.
Can you please tell me how to solve the problem? my goal is: when I give a strain, the RVE returns a correct stress (i.e., the one in the stressStrain.txt )

I would greatly appreciate it if you can help me!
Thank you!