ossillation_rotation

Asked by nobody on 2020-02-14

hi all,

I want to achieve the oscillation movement, as I mentioned in a previous question.[1]

generally speaking, making the shaft rotates 200 degrees in one direction and then reverse it back.

Here is my code.

########################
from yade import ymport
from yade import pack,export,geom
import itertools
from numpy import *
import numpy as np
from yade import plot, export, utils
import math
target=0.12
young=4e8
finalcompFricDegree=19.5
rate = 20*(2*pi/60)
############################################# sphere particles material properties ########################
sphereMat=O.materials.append(CohFrictMat(young=young,poisson=0.3,frictionAngle=radians(finalcompFricDegree),isCohesive=False,alphaKr=0.2,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True,density=2648))
#O.bodies.append(ymport.text("rotation.txt",material = sphereMat))
O.bodies.append([
 sphere((0.1,0.1,0.1),.03,color=(0,1,0)),
])
################################################ create the penetrator and shaft and the cylinder ##############################
facetMat=O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(finalcompFricDegree),density=0))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

facets = O.bodies.append(geom.facetCylinder((200e-3,200e-3,300e-3),200e-3,600e-3,wallMask=6,material='walls',segmentsNumber=100))

x0=0.2;y0=0.2;z0=0.39665
cylinderIDS= O.bodies.append(geom.facetCylinder((x0,y0,z0-0.02),
 radius=0.0125,height=0.15,orientation=Quaternion((0, 0, 1), 0),wallMask=5,segmentsNumber=10, angleRange=None,material=facetMat))
coneIDS= O.bodies.append(geom.facetCone((x0,y0,z0-0.085825-0.02),
 radiusTop=0.0125,radiusBottom=0.0,height=0.02165,orientation=Quaternion((0, 0, 1), 0),wallMask=6,segmentsNumber=10, angleRange=None,material=facetMat))
conecheck=O.bodies[-1].id
################################################ set the global engine for the simulation##############################
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_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),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.7),
 HelixEngine(angularVelocity=rate,linearVelocity=-0.05,rotateAroundZero=True,zeroPoint=(x0,y0,z0),rotationAxis=(0,0,1),ids=cylinderIDS),
 HelixEngine(angularVelocity=rate,linearVelocity=-0.05,rotateAroundZero=True,zeroPoint=(x0,y0,z0),rotationAxis=(0,0,1),ids=coneIDS),
 PyRunner(command='stop_loading()',iterPeriod=10000),
 PyRunner (command ="ossillation_rotation",iterPeriod = 1),
]

def ossillation_rotation():
 rotationAngle = O.bodies[cylinderIDS].state.rot().norm()
 if 0 < rotationAngle <= 200:
  angularVelocity = rate
  O.bodies[cylinderIDS].state.refOri=O.bodies[cylinderIDS].state.ori
 if 0 < rotationAngle <= 200:
  angularVelocity = - rate ### it seems code here is not working.
  O.bodies[cylinderIDS].state.refOri=O.bodies[cylinderIDS].state.ori
########################## stop the simulation ###################################################################
def stop_loading():
 if ((O.bodies[conecheck].state.pos[2]-O.bodies[conecheck].state.refPos[2])*(-1))>target:
  O.pause()

###############################
i defined a function called ossillation_rotation and use pyrunner to check it every timestep. but it seems it doesn'e work. it rotates always in one direction.

Thanks in advance.

Yong

[1] https://answers.launchpad.net/yade/+question/688324

Question information

Language:
English Edit question
Status:
Expired
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-02-20
Last reply:
2020-03-06
Robert Caulk (rcaulk) said : #1

Hello,

You are close, but you want to adjust the member angularVelocity, instead here you are just creating a python variable called angularVelocity.

So add a label to helix

HelixEngine(angularVelocity=rate,linearVelocity=-0.05,rotateAroundZero=True,zeroPoint=(x0,y0,z0),rotationAxis=(0,0,1),ids=coneIDS,label=helix),

then when you want to change the member angularVelocity, call it as:

helix.angularVelocity = newRate

nobody (nobody01) said : #2

Hi Robert,

I did as you mentioned.

here is the code:

############
from yade import ymport
from yade import pack,export,geom
import itertools
from numpy import *
import numpy as np
from yade import plot, export, utils
import math
target=0.12
young=4e8
finalcompFricDegree=19.5
rate = 10*(2*pi/60)
############################################# sphere particles material properties ########################
sphereMat=O.materials.append(CohFrictMat(young=young,poisson=0.3,frictionAngle=radians(finalcompFricDegree),isCohesive=False,alphaKr=0.2,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True,density=2648))
#O.bodies.append(ymport.text("rotation.txt",material = sphereMat))
O.bodies.append([
 sphere((0.1,0.1,0.1),.03,color=(0,1,0)),
])
################################################ create the penetrator and shaft and the cylinder ##############################
facetMat=O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(finalcompFricDegree),density=0))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

facets = O.bodies.append(geom.facetCylinder((200e-3,200e-3,300e-3),200e-3,600e-3,wallMask=6,material='walls',segmentsNumber=100))

x0=0.2;y0=0.2;z0=0.39665
cylinderIDS= O.bodies.append(geom.facetCylinder((x0,y0,z0-0.02),
 radius=0.0125,height=0.15,orientation=Quaternion((0, 0, 1), 0),wallMask=5,segmentsNumber=10, angleRange=None,material=facetMat))
coneIDS= O.bodies.append(geom.facetCone((x0,y0,z0-0.085825-0.02),
 radiusTop=0.0125,radiusBottom=0.0,height=0.02165,orientation=Quaternion((0, 0, 1), 0),wallMask=6,segmentsNumber=10, angleRange=None,material=facetMat))
conecheck=O.bodies[-1].id
################################################ set the global engine for the simulation##############################
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_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),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.7),
 HelixEngine(angularVelocity=rate,linearVelocity=-0.05,rotateAroundZero=True,zeroPoint=(x0,y0,z0),rotationAxis=(0,0,1),ids=cylinderIDS,label ="cylinderRot"),
 HelixEngine(angularVelocity=rate,linearVelocity=-0.05,rotateAroundZero=True,zeroPoint=(x0,y0,z0),rotationAxis=(0,0,1),ids=coneIDS,label = "coneRot"),
 PyRunner(command='stop_loading()',iterPeriod=10000),
 PyRunner (command ="ossillation_rotation()",iterPeriod = 1),
]

def ossillation_rotation():
 rotationAngle = O.bodies[cylinderIDS].state.rot().norm()
 #rotationAngle = cylinderRot.angleTurned
 #print cylinderRot.angleTurned
 #if 0 < rotationAngle <= pi:
 if rotationAngle>0 and rotationAngle <=math.pi:
  #angularVelocity = rate
  cylinderRot.angularVelocity=rate
  coneRot.angularVelocity=rate
  #print cylinderRot.angularVelocity,coneRot.angularVelocity
  O.bodies[cylinderIDS].state.refOri=O.bodies[cylinderIDS].state.ori
  O.bodies[coneIDS].state.refOri=O.bodies[coneIDS].state.ori
 #elif pi < rotationAngle <= 2*pi:
 elif rotationAngle>0 and rotationAngle <= math.pi:
  #angularVelocity = rate
  cylinderRot.angularVelocity=-1*rate
  coneRot.angularVelocity=-1*rate
  #rotationAngle += rate*O.time
  #print cylinderRot.angularVelocity, coneRot.angularVelocity
  #print rotationAngle
  O.bodies[cylinderIDS].state.refOri=O.bodies[cylinderIDS].state.ori
  O.bodies[coneIDS].state.refOri=O.bodies[coneIDS].state.ori
########################## stop the simulation ###################################################################
def stop_loading():
 if ((O.bodies[conecheck].state.pos[2]-O.bodies[conecheck].state.refPos[2])*(-1))>target:
  O.pause()
#O.run()
#O.saveTmp()

But it gives me the error:
ArgumentError: Python argument types in
    BodyContainer.__getitem__(BodyContainer, list)
did not match C++ signature:
    __getitem__(pyBodyContainer {lvalue}, int)

it seems there is a problem for this line: rotationAngle = O.bodies[cylinderIDS].state.rot().norm()

Robert Caulk (rcaulk) said : #3

Hello,

I am happy to help you with your new problem. In the future please post new questions into new threads, as we request in [1]#5.

>>it seems there is a problem for this line: rotationAngle = O.bodies[cylinderIDS].state.rot().norm()

cylinderIDS is a list. You cannot index O.bodies with a list. You can only index with a single integer. I think you can solve the rest on your own :-)

Cheers,

Robert

nobody (nobody01) said : #5

Hi Robert,

I changed the code as you mentioned, right now there is no error.

but it seems the code is not working, the shaft rotates always in one direction. which means the ossillation_rotation function doesn't run at all?
#########################
from yade import ymport
from yade import pack,export,geom
import itertools
from numpy import *
import numpy as np
from yade import plot, export, utils
import math
target=0.12
young=4e8
finalcompFricDegree=19.5
rate = 10*(2*pi/60)
############################################# sphere particles material properties ########################
sphereMat=O.materials.append(CohFrictMat(young=young,poisson=0.3,frictionAngle=radians(finalcompFricDegree),isCohesive=False,alphaKr=0.2,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True,density=2648))
#O.bodies.append(ymport.text("rotation.txt",material = sphereMat))
O.bodies.append([
 sphere((0.1,0.1,0.1),.03,color=(0,1,0)),
])
################################################ create the penetrator and shaft and the cylinder ##############################
facetMat=O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(finalcompFricDegree),density=0))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

#facets = O.bodies.append(geom.facetCylinder((200e-3,200e-3,300e-3),200e-3,600e-3,wallMask=6,material='walls',segmentsNumber=100))

x0=0.2;y0=0.2;z0=0.39665
cylinderIDS= O.bodies.append(geom.facetCylinder((x0,y0,z0-0.02),
 radius=0.0125,height=0.15,orientation=Quaternion((0, 0, 1), 0),wallMask=5,segmentsNumber=10, angleRange=None,material=facetMat))
coneIDS= O.bodies.append(geom.facetCone((x0,y0,z0-0.085825-0.02),
 radiusTop=0.0125,radiusBottom=0.0,height=0.02165,orientation=Quaternion((0, 0, 1), 0),wallMask=6,segmentsNumber=10, angleRange=None,material=facetMat))
cylinder =O.bodies[-2].id
conecheck=O.bodies[-1].id
################################################ set the global engine for the simulation##############################
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_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),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.7),
 HelixEngine(angularVelocity=rate,linearVelocity=-0.05,rotateAroundZero=True,zeroPoint=(x0,y0,z0),rotationAxis=(0,0,1),ids=cylinderIDS,label ="cylinderRot"),
 HelixEngine(angularVelocity=rate,linearVelocity=-0.05,rotateAroundZero=True,zeroPoint=(x0,y0,z0),rotationAxis=(0,0,1),ids=coneIDS,label = "coneRot"),
 PyRunner(command='stop_loading()',iterPeriod=10000),
 PyRunner(command ="ossillation_rotation()",iterPeriod = 10),
]

def ossillation_rotation():
 for b in O.bodies:
  if isinstance(b.shape,Facet):
   rotationAngle = b.state.rot().norm()
   if rotationAngle > 0 and rotationAngle <= math.pi:
    cylinderRot.angularVelocity=rate
    coneRot.angularVelocity=rate
    b.state.refOri=b.state.ori
   elif b.state.rot().norm() > 0 and b.state.rot().norm() <= math.pi:
    cylinderRot.angularVelocity=-1*rate
    coneRot.angularVelocity=-1*rate
    b.state.refOri=b.state.ori
########################## stop the simulation ###################################################################
def stop_loading():
 if ((O.bodies[conecheck].state.pos[2]-O.bodies[conecheck].state.refPos[2])*(-1))>target:
  O.pause()
#O.run()
#O.saveTmp()

##################
can you help me to check why this function doesn't work?

thanks,
Yong

Jan Stránský (honzik) said : #6

> rotationAngle = b.state.rot().norm()
> if rotationAngle > 0 and rotationAngle <= math.pi:
> ...
> elif b.state.rot().norm() > 0 and b.state.rot().norm() <= math.pi:

is if and elif condition different?

cheers
Jan

nobody (nobody01) said : #7

Hi Jan,

I'm thinking when the shaft rotates, the reference orientation will also change.

So I use this line: b.state.refOri=b.state.ori , to set the reference orientation of the end of the timestep.

or we don't need to change the reference orientation at all?

I'm also a little bit confused.

Jan Stránský (honzik) said : #8

my point was, that inside ossillation_rotation you have:
###
if condition:
   cylinderRot.angularVelocity=rate
elif theSameCondition:
   cylinderRot.angularVelocity=-1*rate
###

the problem is that the code inside elif branch is never executed (either it is not fullfiled or is "eaten" by the first if) and the rate is never switched, resulting in "the shaft rotates always in one direction"

cheers
Jan

nobody (nobody01) said : #9

Hi Jan,

sorry for the late reply!

I change the code as you mentioned.
But it still doesn't work. the shaft still rotates in one direction. I think this is related to the reference orientation.
under the elif condition.

if condition:
   cylinderRot.angularVelocity=rate
elif theSameCondition:
   cylinderRot.angularVelocity=-1*rate
If the reference orientation didn't change, for example, when the shaft rotates 200 degrees, if the elif condition work, the shaft will rotate in another direction, maybe with one timestep, the rotation angle will again be 200 degrees. so, we can't see the shaft rotates in another direction.

best,
Yong

Jan Stránský (honzik) said : #10

> I change the code as you mentioned.

please provide the code

cheers
Jan

nobody (nobody01) said : #11

hi Jan,

here is the code:
###########################
from yade import ymport
from yade import pack,export,geom
import itertools
from numpy import *
import numpy as np
from yade import plot, export, utils
import math
target=0.12
young=4e8
finalcompFricDegree=19.5
rate = 10*(2*pi/60)
############################################# sphere particles material properties ########################
sphereMat=O.materials.append(CohFrictMat(young=young,poisson=0.3,frictionAngle=radians(finalcompFricDegree),isCohesive=False,alphaKr=0.2,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True,density=2648))
#O.bodies.append(ymport.text("rotation.txt",material = sphereMat))
O.bodies.append([
 sphere((0.1,0.1,0.1),.03,color=(0,1,0)),
])
################################################ create the penetrator and shaft and the cylinder ##############################
facetMat=O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(finalcompFricDegree),density=0))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

#facets = O.bodies.append(geom.facetCylinder((200e-3,200e-3,300e-3),200e-3,600e-3,wallMask=6,material='walls',segmentsNumber=100))

x0=0.2;y0=0.2;z0=0.39665
cylinderIDS= O.bodies.append(geom.facetCylinder((x0,y0,z0-0.02),
 radius=0.0125,height=0.15,orientation=Quaternion((0, 0, 1), 0),wallMask=5,segmentsNumber=10, angleRange=None,material=facetMat))
coneIDS= O.bodies.append(geom.facetCone((x0,y0,z0-0.085825-0.02),
 radiusTop=0.0125,radiusBottom=0.0,height=0.02165,orientation=Quaternion((0, 0, 1), 0),wallMask=6,segmentsNumber=10, angleRange=None,material=facetMat))
cylinder =O.bodies[-2].id
conecheck=O.bodies[-1].id
################################################ set the global engine for the simulation##############################
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_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),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.7),
 HelixEngine(angularVelocity=rate,linearVelocity=-0.05,rotateAroundZero=True,zeroPoint=(x0,y0,z0),rotationAxis=(0,0,1),ids=cylinderIDS,label ="cylinderRot"),
 HelixEngine(angularVelocity=rate,linearVelocity=-0.05,rotateAroundZero=True,zeroPoint=(x0,y0,z0),rotationAxis=(0,0,1),ids=coneIDS,label = "coneRot"),
 PyRunner(command='stop_loading()',iterPeriod=10000),
 PyRunner(command ="ossillation_rotation()",iterPeriod = 10),
]

def ossillation_rotation():
 for b in O.bodies:
  if isinstance(b.shape,Facet):
   rotationAngle = b.state.rot().norm()
   if rotationAngle > 0 and rotationAngle <= math.pi:
    cylinderRot.angularVelocity=rate
    coneRot.angularVelocity=rate
   elif rotationAngle > 0 and rotationAngle <= math.pi:
    cylinderRot.angularVelocity=-1*rate
    coneRot.angularVelocity=-1*rate
########################## stop the simulation ###################################################################
def stop_loading():
 if ((O.bodies[conecheck].state.pos[2]-O.bodies[conecheck].state.refPos[2])*(-1))>target:
  O.pause()
#O.run()
#O.saveTmp()
####################################
it still rotates in one direction.

best,
Yong

Jan Stránský (honzik) said : #12

> if rotationAngle > 0 and rotationAngle <= math.pi:
> elif rotationAngle > 0 and rotationAngle <= math.pi:

still the same, you have an if condition and exactly the same elif condition. It does not make sense.
If the condition is met, the first if branch is executed and and nothing else.
If the condition is not met, nor it is in elif.
Result = the code in elif branch (setting -1*rate) is NEVER executed -> the rate is never switched, resulting in "the shaft still rotates in one direction".

why not simply else instead of elif? is there any other option?

cheers
Jan

nobody (nobody01) said : #13

Hi Jan,

as you mentioned, I just change the elif to else:

###########
def ossillation_rotation():
 for b in O.bodies:
  if isinstance(b.shape,Facet):
   rotationAngle = b.state.rot().norm()
   if rotationAngle > 0 and rotationAngle <= math.pi:
    cylinderRot.angularVelocity=rate
    coneRot.angularVelocity=rate
   else:
    cylinderRot.angularVelocity=-1*rate
    coneRot.angularVelocity=-1*rate

############
even if I use this, the shaft still rotates in one direction.

best,
Yong

Jan Stránský (honzik) said : #14

Now when the if-elif structure is becoming meaningful, we can continue :-)

with rotation there is a problem that the rotation vector representing the rotation from initial orientation just takes the two states, no history. So if you rotate and angle < pi, it is the rot().
However, if true rotation is > pi, the rot() represents "closest path" and becomes angle-pi and the rot().norm() is never greater than pi (you never comes to the else branch)
###
b = sphere((0,0,0),1)
for f in (0.1,0.5,0.9,1.1,1.5,1.9):
 b.state.ori = Quaternion((1,0,0),f*pi)
 print(f,b.state.rot(),b.state.rot().norm())
###

A solution is to store the rotation incrementally.

But before continuing, why do you have "for b in O.bodies" in ossillation_rotation? Isn't the change of rotation dependent on one specific body?

cheers
Jan

nobody (nobody01) said : #15

Hi Jan,
Because in my simulation I still have some particles, the particles also have the orientations, I just want the shaft rotates under certain degrees. so I use, "for b in O.bodies" and "if isinstance(b.shape,Facet):" to isolate the degrees that the shaft has been rotated.

Today, I also found that the rotation angle will smaller than pi by using O.bodies,state.rot().norm().

so if I want the shaft rotates 200 degrees, I can't use,O.bodies.state.rot().nrom(), this command??

best,
Yong

Launchpad Janitor (janitor) said : #16

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Jan Stránský (honzik) said : #17

Hi,

sorry for late answer..

> Because in my simulation I still have some particles, the particles also have the orientations, I just want the shaft rotates under certain degrees. so I use, "for b in O.bodies" and "if isinstance(b.shape,Facet):" to isolate the degrees that the shaft has been rotated.

in principle it should work as it is now, but the for loop is no necessary and overall it more error prone IMO

cylinderRot and coneRot controls all of its bodies, so it should be sufficient to only use e.g. its first body:
###
def ossillation_rotation():
 i = cylinderRot.ids[0]
 b = O.bodies[i]
 rotationAngle = b.state.rot().norm() # needs to update, but the point is using just one body
 if rotationAngle > 0 ...
  ...
 else:
  ...
###

> so if I want the shaft rotates 200 degrees, I can't use,O.bodies.state.rot().nrom(), this command??

yes, you cannot use it. As I have written, the solution is to store actual rotation and add incrementally the rotation increment to the stored value. MWE with one facet:
###
b = facet(((0,0,0),(1,0,0),(0,1,0)))
O.bodies.append(b)
b.state.angVel = (0,0,.1)

rot = 0
ori = b.state.refOri

def updateRot(b,rot,ori): # works for rotation aroud z axis
   relRot = ori.conjugate()*b.state.ori # same as State.rot(), https://gitlab.com/yade-dev/trunk/-/blob/master/core/State.hpp#L40
   axis,angle = relRot.toAxisAngle() # get axis and angle from quaternion
   if (axis[2]) < 0: # "negative" rotation ...
      angle *= -1 # ... negate the positive angle
   rot += angle # updates rotation angle
   ori = b.state.ori # "store" new orientation value
   return rot,ori

for i in range(20):
   O.step()
   rot,ori = updateRot(b,rot,ori)
   print i,rot

b.state.angVel = (0,0,-.3)

for i in range(40):
   O.step()
   rot,ori = updateRot(b,rot,ori)
   print i,rot
###

cheers
Jan

[1]