How can I control the variables in the Partial engine during the simulation?

Asked by JINSUN LEE

I want to move the facet in a sinusoidal motion.
So, I used the Partial engine(Rotation Engine) and changed angular velocity with virtual time.

For example, O.engines is defined as follows;

O.engines=O.engines+[RotationEngine(rotateAroundZero=True, zeroPoint=(0,0,0), rotationAxis=(1, 0, 0), angularVelocity=rad_vel, ids=- -)]

And the angularVelocity is changed by function, cal_wall_rot()

def cal_wall_rot():
       global rad_vel
       rad_vel = "function of sine(O.time)"

Finally, the function, "cal_wall_rot()" is called in every step.
O.engines=O.engines+[PyRunner(command='cal_wall_rot()', iterPeriod=1)]

However, the variable "rad_vel" is not changed during the simulation.

Anyway, I have found another way to change the velocity of the facet using b.state.

Is this the only way to control the motion of bodies with time?

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
JINSUN LEE
Solved:
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said :
#1

Hello,

First of all, always try to provide a MWE [1], such that we know how exactly the code is used and in which context.

> O.engines=O.engines+[...]

**My** recommendation is not to use this this king of structure unless necessary.
Better approach is to set everything while creating O.engines, i.e.
###
O.engines = [
... # all engines here
]
###

> However, the variable "rad_vel" is not changed during the simulation.

Surely the variable "rad_vel" is changed, just print it:
###
def cal_wall_rot():
       global rad_vel
       rad_vel = "function of sine(O.time)"
       print(rad_vel) # should print current value
###

That it has no effect on the RotationEngine is another point.
rad_vel is a value. If you pass it to RotationEngine(...,angularVelocity=rad_vel), and then you change rad_vel, it has no effect on angularVelocity, because the angularVelocity value is just copy of passed rad_vel at the time of creation.
You need to modify the attribute of the RotationEngine instance, not a global variable rad_vel.

> Anyway, I have found another way to change the velocity of the facet using b.state.
> Is this the only way to control the motion of bodies with time?

yes and no (depending on definition of "the only way")
Modifying directly state is one option.
Using RotationEngine (or other similar engine, like HelixEngine) is another option. However, internally they just modify state.. (hence the dependency on definition)

===============
Solution:
modify attribute of RotationEngine instance, not an independent rad_vel variable.

First, you need to access to the instance. A few options:

- using instance directly
###
rot = RotationEngine(...)
O.engines = [
   ...,
   rot,
]
###

- using label
###
O.engines = [
    ...,
    RotationEngine(...,label="rot"),
]

- using O.engines[index]:
###
rot = O.engines[indexOfRotationEngine]
###

Once having "rot" instance, you can do whatever you like with it, in your case most likely something like:
###
def cal_wall_rot():
       rot.angularVelocity = "function of sine(O.time)"
###

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

Revision history for this message
JINSUN LEE (ppk21) said :
#2

Thank you for the answer. It helps me a lot.

I have two another question,

1. What is the meaning of MWE [1] ?

2. It seems to be a bug when using the label in 'RotationEngine'

###
In [1]: O.engines
Out[1]:
[<ForceResetter instance at 0x559595ddcf40>,
 <InsertionSortCollider instance at 0x559595f91560>,
 <InteractionLoop instance at 0x559596078500>,
 <NewtonIntegrator instance at 0x559595f9d0c0>,
 <PyRunner instance at 0x559595c04200>,
 <PyRunner instance at 0x559595b067a0>,
 <ServoPIDController instance at 0x559595a9ac20>,
 <DomainLimiter instance at 0x559595f195c0>,
 <PyRunner instance at 0x559595b06ba0>,
 <RotationEngine instance at 0x559596019160>]
###

When I check the label of the partial engine, label can be identified
The label of ServoPIDController is instance,
However, the label of RontationEngine gives error message as follows,

###
In [2]: O.engines[6].label
Out[2]: 'servo_1'

In [3]: servo_1
Out[3]: <ServoPIDController instance at 0x559595a9ac20>

In [4]: O.engines[9].label
Out[4]: 'moving_wall_front'

In [5]: moving_wall_front
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
/usr/bin/yade in <module>
----> 1 moving_wall_front

NameError: name 'moving_wall_front' is not defined
###

Ubuntu 22.10

Welcome to Yade 2022.01a
Using python version: 3.10.7 (main, Nov 24 2022, 19:45:47) [GCC 12.2.0]

Revision history for this message
Jan Stránský (honzik) said :
#3

> 1. What is the meaning of MWE [1] ?

Meaning of MWE is explained in [1] :-)
[1] means a link, for convenience not put directly in the text, but referenced by number as "[1]" and listed at the end of the message, e.g. as
[1] https://www.yade-dem.org/wiki/Howtoask
So in this case you should go to https://www.yade-dem.org/wiki/Howtoask and read about MWE.

> 2. It seems to be a bug when using the label in 'RotationEngine'

please provide a MWE [1], without it we can merely guess what is going on..

Cheers
Jan

Revision history for this message
JINSUN LEE (ppk21) said :
#4

 Dear, Jan.

 Thank you for your kind and prompt reply.

 I tried to make MWE as simple as it is. but it contains local PATH.
 So, please check the line "sys.path.append('/home/jinsun/Dropbox/Yadee/Test')"

 The script consist of 2 sequential script (c.py, c2.py)

You can find the follwoing error after running the final script (c2.py)

The operating system is UBUNTU 22.10

===
% yade c2.py
Welcome to Yade 2022.01a
Using python version: 3.10.7 (main, Nov 24 2022, 19:45:47) [GCC 12.2.0]
TCP python prompt on localhost:9000, auth cookie `kdscae'
XMLRPC info provider on http://localhost:21000
Running script c2.py
[[ ^L clears screen, ^U kills line. F12 controller, F11 3D view (press "h" in 3D view for help), F10 both, F9 generator, F8 plot. ]]

In [1]: O.engines
Out[1]:
[<ForceResetter instance at 0x55a9d9360d90>,
 <InsertionSortCollider instance at 0x55a9d943a220>,
 <InteractionLoop instance at 0x55a9d930dff0>,
 <NewtonIntegrator instance at 0x55a9d9120d20>,
 <ServoPIDController instance at 0x55a9d95575a0>,
 <PyRunner instance at 0x55a9d90a7dc0>,
 <RotationEngine instance at 0x55a9d91173d0>]

In [2]: O.engines[6].label
Out[2]: 'ak'

In [3]: ak.ids
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
/usr/bin/yade in <module>
----> 1 ak.ids

NameError: name 'ak' is not defined
===

Script c.py
===
from yade import pack, plot

def setGeomVars (): # initialize variables
     width = 0.05
     height = 0.03
     margin = 30
     # Calculate extra length
     dx = width/100/2*margin
     dy = width/100/2*margin
     dz = height/100/2*margin
     saveVars('geoms',loadnow=True,**locals())

setGeomVars ()
from yade.params.geoms import * # load initilized variables

# side pannel
p1s = (-width/2,-width/2-dy,-height/2-dz)
p5s = (-width/2,-width/2-dy,height/2+dz)
p6s = (-width/2,width/2+dy,height/2+dz)
p2s = (-width/2,width/2+dy,-height/2-dz)

side1_1 = utils.facet(vertices=[p1s,p5s,p2s], wire=True, highlight=False) #p1 p5 p2
side1_2 = utils.facet(vertices=[p2s,p5s,p6s], wire=True, highlight=False) #p1 p5 p2
O.bodies.append(side1_1)
O.bodies.append(side1_2)

p4s = (width/2,-width/2-dy,-height/2-dz)
p8s = (width/2,-width/2-dy,height/2+dz)
p7s = (width/2,width/2+dy,height/2+dz)
p3s = (width/2,width/2+dy,-height/2-dz)

side2_1 = utils.facet(vertices=[p4s,p8s,p3s], wire=True, highlight=False) #p1 p5 p2
side2_2 = utils.facet(vertices=[p3s,p8s,p7s], wire=True, highlight=False) #p1 p5 p2
O.bodies.append(side2_1)
O.bodies.append(side2_2)

cp1 = (-width/2,-width/2,2*height/2+dz) #corner point #1
cp2 = (width/2,width/2,3*height/2+dz) #corner point #2
radius_mean = 0.001

sp = pack.SpherePack()
sp.makeCloud(cp1, cp2, rMean=radius_mean, rRelFuzz=0.0, num = 1000)
sp.toSimulation(color=(0,0,1)) # pure blue

O.materials.append(FrictMat(young=20e6, poisson=0.17, density=2000))

wallid1 = range(0,2)

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
    InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()], # collision geometry
            [Ip2_FrictMat_FrictMat_FrictPhys()], # collision "physics"
            [Law2_ScGeom_FrictPhys_CundallStrack()] # contact law -- apply forces
    ),
    NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.1),
    ServoPIDController(axis=(0,0,1), ids=wallid1, target=1000, kD=5.0, kI=5.0, kP=5.0, maxVelocity=-0.01, iterPeriod=100, label='ab')
]

def addData():
     for b in O.bodies:
         b.shape.color=scalarOnColorScale(b.state.vel.norm())

O.run(100)

O.engines=O.engines+[PyRunner(command='addData()', iterPeriod=100)]

O.run(100)

O.save("c1.yade")
===

Script c2.py
===
from yade import pack, plot
import sys
import math
sys.path.append('/home/jinsun/Dropbox/Yadee/Test')
O.load("c1.yade")
loadVars('geoms')

from yade.params.geoms import * # load initilized variables

wallid2 = range(2,4)

def addData():
     for b in O.bodies:
         b.shape.color=scalarOnColorScale(b.state.vel.norm())

O.engines=O.engines+[RotationEngine(rotateAroundZero=True, zeroPoint=(0,0,0), rotationAxis=(1, 0, 0), angularVelocity=0.0, ids=wallid2, label="ak")]
===

END OF THE SCRIPT

Revision history for this message
Jan Stránský (honzik) said :
#5

> O.run(100)
> O.engines=O.engines+[PyRunner(command='addData()', iterPeriod=100)]
> O.run(100)
> O.save("c1.yade")

Use O.run(100,True), or more verbose O.run(100,wait=True). See [2]. Then the run command wait until 100 steps are finished (probably meant like that).
If I used it this way, the labels worked.
I have no idea why without wait=True it does not work..

You can also use the other approaches mentioned in my first answer:
- instance appraoch: ak=RotationEngine(...)
- index approach: ak=O.engines[index]

Cheers
Jan

[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Omega.run

Revision history for this message
JINSUN LEE (ppk21) said :
#6

Thx. Appreciate :)