forcerecorder and O.force.t does not give correct value

Hi,

I am trying to insert a big clump into a bed of spheres, displace this clump in the spheres and record force and torque on this clump.
However, the torque I get does not seem to make sense.

I tried both O.forces.t and torque recorder.
The O.forces.t gives torque that does not match the force on the clump (the force is about 100 N and 0.15 meter away from the center of the clump but the torque recorded is only about 5N*m.

The torque recorder gives values that does not make much sense. I guess the value it recorded depends on the rotation axis and zeroPoint. How are these defined w.r.t to? Are they in the global coordinates or the clump's local coordinates?

Below is the code if you want to take a look at it.

Thanks,
Yifan

#!/usr/bin/python
# -*- coding: utf-8 -*-

import random
from yade import ymport
from yade import qt
from yade import plot
import sys,time

factor=1
test='2'
level=4700*factor
## PhysicalParameters
sand=O.materials.append(FrictMat(young=3e8,poisson=0.3,density=1631,frictionAngle=0.5236))
#default=O.materials.append(FrictMat(young=1e7,poisson=.3,density=1e3,frictionAngle=.5))
default=O.materials.append(CohFrictMat(young=1e7,poisson=.3,density=1e3,frictionAngle=.5,alphaKr=2.0,alphaKtw=1.0,etaRoll=-1,isCohesive=False,momentRotationLaw=True))
# Clump
clump1=O.bodies.appendClumped([\
sphere([0,0,0.1], radius=0.1,material=default),\
sphere([0,0,0.15],radius=0.1,material=default),\
sphere([0,0,0.2], radius=0.1,material=default),\
sphere([0,0,0.25], radius=0.1,material=default),\
sphere([0,0,0.3], radius=0.1,material=default),\
sphere([0,0,0.35], radius=0.1,material=default),\
sphere([0,0,0.4], radius=0.1,material=default),\
])

#s=O.bodies.append(sphere(center=(0,0,0.5),radius=0.2,fixed=False))
#get clump ids:
ClumpID = clump1[0]
SpheresID = clump1[1]
O.bodies[ClumpID].state.blockedDOFs='XYZ'
O.bodies[ClumpID].state.angVel=(0,0,0)
O.bodies[ClumpID].state.vel=(0,0,0)
#print id_clump2
#Default Material:
#FrictMat(density=1e3,young=1e7,poisson=.3,frictionAngle=.5).

# Spheres
sphereRadius = 0.01
#nbSpheres = (45,45,60)
nbSpheres=(32,32,30)
# with 32,32,11 spheres settle around -0.08555m after 2s

surfacesphere=[]
print "Creating %d spheres..."%(nbSpheres[0]*nbSpheres[1]*nbSpheres[2]),
for i in xrange(nbSpheres[0]):
    for j in xrange(nbSpheres[1]):
  for k in xrange(nbSpheres[2]):
   x = (i*2 - nbSpheres[0])*sphereRadius*1.1+sphereRadius*random.uniform(-0.1,0.1)
   y = (j*2 - nbSpheres[1])*sphereRadius*1.1+sphereRadius*random.uniform(-0.1,0.1)
   z = -k*sphereRadius*2.2-0.01
   r = random.uniform(sphereRadius,sphereRadius*0.9)
   fixed = False
   color=[0.51,0.52,0.4]
   if (i==0 or i==nbSpheres[0]-1 or j==nbSpheres[1]-1 or j==0 or k==nbSpheres[2]-1):
    fixed = True
    color=[0.21,0.22,0.1]
   tmp=O.bodies.append(sphere([x,y,z],r,color=color,fixed=fixed,material=default))
   if k==0 and not fixed:
    surfacesphere.append(tmp)

print "done\n"

## Estimate time step
#O.dt=PWaveTimeStep()
O.dt=0.0001

def changeforce():
 if O.iter <15000*factor:
  O.bodies[ClumpID].state.pos=(-0.1,0,0.25)
  O.bodies[ClumpID].state.angVel=(0,0,0)
  O.bodies[ClumpID].state.vel=(0,0,0)
 elif O.iter == 15000*factor:
  O.bodies[ClumpID].state.pos=(-0.1,0,0.25-0.238)
  O.bodies[ClumpID].state.angVel=(0,0,0)
  O.bodies[ClumpID].state.vel=(0,0,0)
 elif O.iter > 15000*factor and O.iter < 15000*factor+16000:
  O.bodies[ClumpID].state.pos=(-0.1,0,0.25-0.238-O.dt*0.1*(O.iter-15000))
 elif O.iter >= 15000*factor+16000:
  O.bodies[ClumpID].state.pos=(-0.1+0.1414*O.dt*(O.iter-31000),0,0.25-0.238-0.16)

def addPlotdata():
 plot.addData(i=O.iter,x=O.bodies[ClumpID].state.pos[0],z=O.bodies[ClumpID].state.pos[2],f=O.forces.f(ClumpID,sync=True),t=O.forces.t(ClumpID,sync=True))
 #print O.bodies[ClumpID].state.pos[2]

def savedata():
 #pass
 if O.iter > 40000*factor:
  plot.saveDataTxt('displacement'+test+'.dat')

O.engines=[
 ## Resets forces and momenta the act on bodies
 ForceResetter(),

 ##display all kinds of info here

 ## Using bounding boxes find possible body collisions.
 InsertionSortCollider([
  Bo1_Sphere_Aabb(),
  Bo1_Facet_Aabb(),
 ]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
  #[Ip2_FrictMat_FrictMat_FrictPhys()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 NewtonIntegrator(damping=0.3,gravity=[0,0,-9.81]),

 ## Save force onClump
 ForceRecorder(ids=[ClumpID],file='force'+test+'.dat',iterPeriod=1),
 TorqueRecorder(ids=[ClumpID],file='itorque'+test+'.dat',iterPeriod=1,rotationAxis=(0,1,0),zeroPoint=(0,0,1)),
 TorqueRecorder(ids=[ClumpID],file='iitorque'+test+'.dat',iterPeriod=1,rotationAxis=(0,1,0),zeroPoint=(0,0,0)),
 TorqueRecorder(ids=[ClumpID],file='iiitorque'+test+'.dat',iterPeriod=1,rotationAxis=(0,1,0),zeroPoint=(0,0,-1)),
 PyRunner(command='changeforce()',iterPeriod = 1),
 PyRunner(command='addPlotdata()',iterPeriod = 1),
 PyRunner(command='savedata()',iterPeriod = 100)
]

nbIter=50000*factor

qt.View()

O.stopAtIter=nbIter
O.run()