How to export relative velocity of spheres to Paraview

Hello Chien-Cheng,

> the relative velocity of all interacted spheres.

relative velocity with respect to what? relative velocity of interacting bodies for all interactions?

you can use export.VTKExporter [1]
###
vtk = export.VTKExporter()
vtk.exportSpheres(what=dict(velRelative="b.state.vel"))
vtk.exportInteractions(what=dict(velRelative="i.phys.normalForce"))
###
above, b.state.vel and i.phys.normalForce is only illustrative, needs to be replaced by actual formula for relative velocity
If you use older Yade version, the correct format might be list of tuples instead of dict: what=[("velRelative","b.state.vel")]

cheers
Jan

[1] https://yade-dem.org/doc/yade.export.html#yade.export.VTKExporter

Hi Jan,

Thanks for your quick reply!

> relative velocity with respect to what? relative velocity of interacting bodies for all interactions?

Sorry I didn't describe clearly. But, yes, the relative velocity of interacting bodies for all interactions.
The relative velocity of the spheres that I want to obtain is described in ScGeom [1].

> above, b.state.vel and i.phys.normalForce is only illustrative, needs to be replaced by actual formula for relative velocity.

Can I get the formula for relative velocity [1] somewhere in Yade? Or I have to define one by myself?

Cheers,
Chien-Cheng

[1] https://yade-dem.org/doc/yade.wrapper.html?highlight=scgeom#yade.wrapper.ScGeom

> Can I get the formula for relative velocity [1] somewhere in Yade? Or I have to define one by myself?

AFAIK you have to write it yourself..

you can write a function to it and call the function in exportInteractions / exportContactPoints, but as the commands are evaluated only locally, you have to "make the function global", e.g. using __builtin__ module (be very careful with this workaround not to overwrite some python basics, e.g. list etc.):
###
def relativeVelocity(interaction):
   i1,i2 = interaction.id1, interaction.id2
   b1,b2 = [O.bodies[i] for i in (i1,i2)]
   r1,r2 = [b.shape.radius for b in (b1,b2)]
   v1,v2 = [b.state.vel for b in (b1,b2)]
   ... # actual code comes here
   return relativeVelocity

import __builtin__ # or buitlins (2ithout underscores) in Python3
__builtin__.relativeVelocity = relativeVelocity

vtk.exportInteractions(what=..."relativeVelocity(i)"...) # now you can use relativeVelocity inside exportInteractions
###

cheers
Jan

Hi Jan,

Thanks a lot!

> you can write a function to it and call the function in exportInteractions / exportContactPoints

Is it possible that instead of calling the function, I can make a list of data of relative velocity that is recognized in Paraview

Cheers,
Chien-Cheng

> Is it possible that instead of calling the function, I can make a list of data of relative velocity that is recognized in Paraview

I am not sure i I got it correctly, but you can e.g. create a dict with (id1,id2) keys and read from it in exportInteractions, something like (not tested):
###
interactionsInfo = {}
for i in O.interactions:
   ... # some computations
   interactionsInfo[(i.id1,i.id2)] = someData
...
vtk.exportInteractions(what=dict(something="interactionsInfo[(i.id1,i.id2)]"))
# probably interactionsInfo would need the same builtin trick
###

cheers
Jan

Thanks Jan Stránský, that solved my question.

Hi Jan,

I would like to export the function in vtk files as below. I adjust the body of function unlike the above you showed me.
But it shows me that "name i is not defined". How do I modify my code?
Thank you for your help again!

###
def weakeningTemp(interaction):
    heatCapacity = 900 # [J/kg*K]
    density = 2500 # [kg/m3]
    totalHeat = 0
    thermal_diff = 1.e-6
    shear_modulus = 22e9
    tau_c = 0.1*shear_modulus
    initialTemp = 20
    for i in O.interactions():
        if O.bodies[i.id1].isReal and O.bodies[i.id2].isReal:
            penetrationDepth = i.geom.penetrationDepth
            radius1 = i.geom.refR1
            radius2 = i.geom.refR2
            effectiveRadius = (i.geom.refR1 * i.geom.refR2 )/( i.geom.refR1+i.geom.refR2 )
            contactRadius = math.sqrt( i.geom.penetrationDepth *effectiveRadius)
            Da = 2 * contactRadius
            contactArea = np.pi*(contactRadius**2)
            relativeVelocity = i.geom.shearInc.norm()/O.dt
            weakeningTemp = math.sqrt((relativeVelocity*Da)/(np.pi*thermal_diff))*tau_c/(density*heatCapacity) + initialTemp
    return weakeningTemp

builtins.weakeningTemp = weakeningTemp # such that VTKExporter has access to it

vtk = export.VTKExporter("weakingTemp")
def vtkExport():
    vtk.exportSpheres(ids=list(range(2,len(O.bodies))),what=dict(weakeningTemp="weakeningTemp(i)"))
###

Cheers,
Chien-Cheng

> But it shows me that "name i is not defined".

please post the complete error message

thanks
Jan

> what=dict(weakeningTemp="weakeningTemp(i)")

weakeningTemp(i) means it is a function of **one interaction**. In your function, the interaction argument is never used and instead you do "for i in O.interactions" instead, which should not be like this..
Should be something like:
###
def weakeningTemp(interaction):
    heatCapacity = 900 # [J/kg*K]
    ...
    initialTemp = 20
    penetrationDepth = i.geom.penetrationDepth
    ...
    weakeningTemp = ...
    return weakeningTemp
###

> if O.bodies[i.id1].isReal and O.bodies[i.id2].isReal:

body does not have isReal member..

cheers
Jan

Hi Jan,

Thanks for your reply.

I still get the same error message after following your instruction.

The complete error message:
###
NameError Traceback (most recent call last)
/home/hung/YadeNov2019/bins/yade-2019-11-18.git-b1a2767 in <module>()

/home/hung/YadeNov2019/bins/yade-2019-11-18.git-b1a2767 in vtkExport()
    183 vtk = export.VTKExporter("weakingTemp")
    184 def vtkExport():
--> 185 vtk.exportSpheres(ids=list(range(2,len(O.bodies))),what=dict(weakeningTemp="weakeningTemp(i)"))
    186
    187 ##################################################

/home/hung/YadeNov2019/build/lib/x86_64-linux-gnu/yade-2019-11-18.git-b1a2767/py/yade/export.py in exportSpheres(self, ids, what, comment, numLabel, useRef)
    427 # write additional data from 'what' param
    428 for name,command in what.items(): # for each name...
--> 429 test = eval(command) # ... eval one example to see what type (float, Vector3, Matrix3) the result is ...
    430 # ... and write appropriate header line and loop over all bodies and write appropriate vtk line(s)
    431 if isinstance(test,Matrix3):

/home/hung/YadeNov2019/build/lib/x86_64-linux-gnu/yade-2019-11-18.git-b1a2767/py/yade/export.py in <module>()

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

Part of my code looks like:
###
O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_FrictMat_FrictMat_MindlinPhys()],
  [Law2_ScGeom_MindlinPhys_Mindlin()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=TSSC,defaultDt=utils.PWaveTimeStep())
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ,PyRunner(command='dataRecorder()',iterPeriod=10,label='recData',dead=True)
 ,PyRunner(command='fixVelocity(RATE_shear)',iterPeriod=1,label='fixVel',dead=True)
 ,PyRunner(iterPeriod=1,command="weakeningTemp",label='weakenTemp',dead=True)
 ,PyRunner(iterPeriod=1,command="vtkExport()",label='vtkExp',dead=True)
 ]

def weakeningTemp(interaction):
    heatCapacity = 900 # [J/kg*K]
    density = 2500 # [kg/m3]
    totalHeat = 0
    thermal_diff = 1.e-6
    shear_modulus = 22e9
    tau_c = 0.1*shear_modulus
    initialTemp = 20
    penetrationDepth = i.geom.penetrationDepth
    radius1 = i.geom.refR1
    radius2 = i.geom.refR2
    effectiveRadius = (i.geom.refR1 * i.geom.refR2 )/( i.geom.refR1+i.geom.refR2 )
    contactRadius = math.sqrt( i.geom.penetrationDepth *effectiveRadius)
    Da = 2 * contactRadius
    contactArea = np.pi*(contactRadius**2)
    relativeVelocity = i.geom.shearInc.norm()/O.dt
    weakeningTemp = math.sqrt((relativeVelocity*Da)/(np.pi*thermal_diff))*tau_c/(density*heatCapacity) + initialTemp
    return weakeningTemp

builtins.weakeningTemp = weakeningTemp # such that VTKExporter has access to it

vtk = export.VTKExporter("weakingTemp")
def vtkExport():
    vtk.exportSpheres(ids=list(range(2,len(O.bodies))),what=dict(weakeningTemp="weakeningTemp(i)"))
###

Cheers,
Chien-Cheng

Thanks for more info

> NameError: name 'i' is not defined
> vtk.exportSpheres(ids=list(range(2,len(O.bodies))),what=dict(weakeningTemp="weakeningTemp(i)"))

vtk.exportSpheres expect "b" as the body.
vtk.exportInteractions expect "i" as the interaction.
You have to choose :-)

cheers
Jan

Thanks Jan Stránský, that solved my question.

Hi Jan,

In my "vtk.exportInteractions" command, I don't want to export the interactions that are related to O.bodies[0] and O.bodies[1]. How do exclude specific interactions from my export vtk files?

I tried as follow but it seems not to work.
I'd appreciate it if you could give some hints. Thank you!

###
def weakeningTemp(i):
    heatCapacity = 900 # [J/kg*K]
    density = 2500 # [kg/m3]
    thermal_diff = 1.e-6
    shear_modulus = 22e9
    tau_c = 0.1*shear_modulus
    initialTemp = 20
    penetrationDepth = i.geom.penetrationDepth
    radius1 = i.geom.refR1
    radius2 = i.geom.refR2
    effectiveRadius = (i.geom.refR1 * i.geom.refR2 )/( i.geom.refR1+i.geom.refR2 )
    contactRadius = math.sqrt( i.geom.penetrationDepth *effectiveRadius)
    Da = 2 * contactRadius
    contactArea = np.pi*(contactRadius**2)
    relativeVelocity = i.geom.shearInc.norm()/O.dt
    weakeningTemp = (math.sqrt((relativeVelocity*Da)/(np.pi*thermal_diff))*tau_c)/(density*heatCapacity) + initialTemp
    return weakeningTemp

builtins.weakeningTemp = weakeningTemp # such that VTKExporter has access to it

vtk = export.VTKExporter("weakingTemp_test")
def vtkExport():
    for a in range(2,len(O.bodies)):
        b = O.bodies[a]
    for i in b.intrs():
        O.interactions.erase(i.id1,0)
        O.interactions.erase(i.id1,1)
        O.interactions.erase(0,i.id2)
        O.interactions.erase(1,i.id2)
        vtk.exportInteractions(what=dict(weakeningTemp="weakeningTemp(i)"))
###

Cheers,
Chien-Cheng

paremter 'ids' [1] is for this:
###
ids = [(i.id1,i.id2) for i in O.interactions] # all ids
ids = [(i,j) for i,j in ids if not i in (0,1) and not j in (0,1)] # filtered ids not containing 0 and 1
vtk.exportInteractions(ids=ids,...)
###

cheers
Jan

[1] https://yade-dem.org/doc/yade.export.html#yade.export.VTKExporter.exportInteractions

Thanks, Jan. That's really helpful!