# ossillation_rotation

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.

########################
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 ########################
#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 ##############################
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

x0=0.2;y0=0.2;z0=0.39665
cylinderIDS= O.bodies.append(geom.facetCylinder((x0,y0,z0-0.02),
coneIDS= O.bodies.append(geom.facetCone((x0,y0,z0-0.085825-0.02),
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 ="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 ###################################################################
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.

Yong

## Question information

Language:
English Edit question
Status:
Expired
For:
Assignee:
No assignee Edit question
Last query:
2020-02-20
2020-03-06
 Robert Caulk (rcaulk) said on 2020-02-14: #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 on 2020-02-14: #2

Hi Robert,

I did as you mentioned.

here is the code:

############
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 ########################
#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 ##############################
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

x0=0.2;y0=0.2;z0=0.39665
cylinderIDS= O.bodies.append(geom.facetCylinder((x0,y0,z0-0.02),
coneIDS= O.bodies.append(geom.facetCone((x0,y0,z0-0.085825-0.02),
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 ="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 ###################################################################
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 on 2020-02-14: #3

Hello,

>>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

 Robert Caulk (rcaulk) said on 2020-02-14: #4
 nobody (nobody01) said on 2020-02-14: #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?
#########################
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 ########################
#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 ##############################
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

x0=0.2;y0=0.2;z0=0.39665
cylinderIDS= O.bodies.append(geom.facetCylinder((x0,y0,z0-0.02),
coneIDS= O.bodies.append(geom.facetCone((x0,y0,z0-0.085825-0.02),
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 ="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 ###################################################################
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 on 2020-02-14: #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 on 2020-02-14: #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 on 2020-02-14: #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 on 2020-02-17: #9

Hi Jan,

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 on 2020-02-17: #10

> I change the code as you mentioned.

cheers
Jan

 nobody (nobody01) said on 2020-02-17: #11

hi Jan,

here is the code:
###########################
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 ########################
#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 ##############################
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

x0=0.2;y0=0.2;z0=0.39665
cylinderIDS= O.bodies.append(geom.facetCylinder((x0,y0,z0-0.02),
coneIDS= O.bodies.append(geom.facetCone((x0,y0,z0-0.085825-0.02),
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 ="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 ###################################################################
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 on 2020-02-17: #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 on 2020-02-18: #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 on 2020-02-18: #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 on 2020-02-20: #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 on 2020-03-06: #16

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

 Jan Stránský (honzik) said on 2020-03-06: #17

Hi,

> 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]