Filterpack predicate definition issue

Asked by Yuxuan Wang on 2020-07-13

Hello Yade programmers and users,
I encountered a problem in my recent simulation script which has been troubling me, and it’d be much appreciated if I can get some advice.

Here are the steps I’m trying to accomplish:
1. Make a loose pack of particles following defined PSD.
2. Use compression engine to generate a denser pack.
3. At defined iteration #, turn off compression engine.
4. Filter pack with pack.filterpack into desired shape.
5. Continue with further steps (gravity deposition, geometry rotation, etc).

Problem:
My error came from step 4, in the "gravityDeposition()" portion of my code. When I tried using pack.filterpack with the engine running, the filter was not applied. I then paused O.engines to check if my predicate definition was truly recognized, and it turns out that it wasn’t for some reason. Would anyone be able to help me find the cause of this? Thanks a lot!

Error:
---------------------------------------------------
NameError Traceback (most recent call last)
/usr/bin/yade in <module>()
----> 1 pred

NameError: name 'pred' is not defined
---------------------------------------------------

Here is my script using yade 2018.02b:
-------------------------
#!/usr/bin/python
# -*- coding: utf-8 -*-

from __future__ import print_function
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab
from yade import pack, plot, qt
from yade import export

### DEFINING VARIABLES AND MATERIALS ###

#geometry definition
Ri,Rs,Ro=6.5,8.5,10.5
H=3.5
p=H/2
k=0.01
Hs=H*k
hs=Hs/2
Rx=2*Ro
Hx=2*H

c1=geom.facetCylinder(center=(0,0,p),radius=Ro, height=H, segmentsNumber=20, wallMask=6)
c2=geom.facetCylinder(center=(0,0,hs),radius=Rs, height=Hs, segmentsNumber=20, wallMask=7)
c3=geom.facetCylinder(center=(0,0,p),radius=Ri, height=H, segmentsNumber=20, wallMask=7)

nRead=readParamsFromTable(
 num_spheres=5000,
 compFricDegree = 30,
 unknownOk=True
)

from yade.params import table
num_spheres=table.num_spheres
compFricDegree = table.compFricDegree
damp=0.2
young=5e6 # contact stiffness
mn,mx=Vector3(-50,-50,-50),Vector3(50,50,50) # corners of the initial packing

## create materials for spheres and plates
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a psd loose particles packing
global sp
sp=pack.SpherePack()
sp.makeCloud(mn,mx,num=6000,psdSizes=[1,2,2.25,2.5,3,4,6],psdCumm=[0.,0.1,0.3,0.3,.3,.7,1.])
sp.toSimulation(material='spheres')

### DEFINING ENGINES ###

triax=TriaxialStressController(
 thickness = 0,
 stressMask = 7,
 internalCompaction=False,
 label='controller'
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),

 PyRunner(command='stop()',iterPeriod=10,label='checker'),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

### APPLYING CONFINING PRESSURE ###

triax.goal1=triax.goal2=triax.goal3=-10000

def stop():
 if O.iter<100:return
 checker.command='gravityDeposition()'

def gravityDeposition():
  pred=pack.inCylinder((0,0,0),(0,0,H),Ro)-pack.inCylinder((0,0,0),(0,0,Hs),Rs)-pack.inCylinder((0,0,0),(0,0,H),Ri)
  assembly=pack.filterSpherePack(pred,sp,True)
  assembly.toSimulation()

#...connect to next steps in simulation
---------------------------------------------------

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-08-02
Last query:
2020-08-02
Last reply:
2020-07-16
Jan Stránský (honzik) said : #1

> When I tried using pack.filterpack with the engine running, the filter was not applied. I then paused O.engines to check if my predicate definition was truly recognized, and it turns out that it wasn’t for some reason.

please be more specific how yo know it was not applied / was not truly recognized?
I have tried your code with different observarions:
###
def gravityDeposition():
  O.pause()
  pred=pack.inCylinder((0,0,0),(0,0,H),Ro)-pack.inCylinder((0,0,0),(0,0,Hs),Rs)-pack.inCylinder((0,0,0),(0,0,H),Ri)
  assembly=pack.filterSpherePack(pred,sp,True)
  print(len(O.bodies))
  assembly.toSimulation()
  print(len(O.bodies))
###
I had 6006 bodies, before, 6008 after assembly.toSimulation()

cheers
Jan

Yuxuan Wang (ywang1) said : #2

Thanks for getting back to me, Jan!
When I entered O.pause() in the terminal and called for "pred", the output was as in my original question, "NameError: name 'pred' is not defined". The same issue happened when I called for "assembly". That was why I thought the predicate was not applied.

Another thing that led me to that conclusion was this: particles in the compressed pack (sp) outside of the predicate (pred) were not deleted after I applied "filterpack". Now coming to think of it, filterpack probably only selects particles within the predicate region but doesn't delete other particles?

However, the fact that you were able to see a reduction in # of bodies means "pred" was recognized by the program. Is that truely the case even though the terminal outputs show they are undefined?

Regards,
Yuxuan

Best Jan Stránský (honzik) said : #3

> When I entered O.pause() in the terminal and called for "pred"

if you call also "pred" form terminal, it is not defined. "pred" and "assembly" are local names defined only inside grevityDeposition function.
So if you want to experiment, add code to that function.

> Another thing ...

there are a few issues in your approach:

1) "sp" is still the same. If you compress particles in the simulation, it has no backward effect on the original "sp". You can do something like:
###
sp.toSimulation()
compression()
sp.fromSimulation()
###
to update "sp"

2) filterSpherePack returns a new sphere pack, again with no influence on the original sphere pack nor the influence on actual simulation. So probably you want something like
###
sp = filterSpherePack(pred,sp,True)
O.bodies.clear() # otherwise all existing particles still exist
sp.toSimulation()
###

> the fact that you were able to see a reduction in # of bodies

no, I saw increase of # of bodies, as your code (considering all above) actually should have done :-)

cheers
Jan

Yuxuan Wang (ywang1) said : #4

Hi Jan,
Thank you for the suggestions and they worked quite well for me upon implementation. I saw a significant reduction of particles going from the compressed pack to the filtered pack.

There is only one thing that I'm still confused about now: I was not able to see any of the particles after "assembly.toSimulation()" in the simulation view, even though I could output their locations. Would you mind giving me a hand on that again?

Thanks again!
Best,
Yuxuan

---------------------
#!/usr/bin/python
# -*- coding: utf-8 -*-
#modifications from trial 10: incorporating Jan's advice on filterpack.

from __future__ import print_function
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab
from yade import pack, plot, qt
from yade import export

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

#geometry definition
Ri,Rs,Ro=6.5,8.5,10.5
H=3.5
p=H/2
k=0.01
Hs=H*k
hs=Hs/2
Rx=2*Ro
Hx=2*H

c1=geom.facetCylinder(center=(0,0,p),radius=Ro, height=H, segmentsNumber=20, wallMask=6)
c2=geom.facetCylinder(center=(0,0,hs),radius=Rs, height=Hs, segmentsNumber=20, wallMask=7)
c3=geom.facetCylinder(center=(0,0,p),radius=Ri, height=H, segmentsNumber=20, wallMask=7)

nRead=readParamsFromTable(
 num_spheres=5000,
 compFricDegree = 30,
 unknownOk=True
)

from yade.params import table
num_spheres=table.num_spheres
compFricDegree = table.compFricDegree
damp=0.2
young=5e6 # contact stiffness
mn,mx=Vector3(-12,-12,-12),Vector3(12,12,12) # corners of the initial packing

## create materials for spheres and plates
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a psd loose particles packing
global sp
sp=pack.SpherePack()
sp.makeCloud(mn,mx,num=1000,psdSizes=[0.399,0.4,0.799,0.8],psdCumm=[0,0.5,0.5,1])
sp.toSimulation(material='spheres')

############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
 thickness = 0,
 stressMask = 7,
 internalCompaction=False,
 label='controller'
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),

 PyRunner(command='stop()',iterPeriod=1,label='checker'),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

#######################################
### APPLYING CONFINING PRESSURE ###
#######################################

triax.goal1=triax.goal2=triax.goal3=-100

def stop():
 if O.iter<100:return
 else:
  controller.dead=True
  for wall in walls:
   wall.state.vel = Vector3.Zero
  for b in O.bodies:
   b.state.vel = Vector3.Zero
  sp.fromSimulation() # update sp from original loose pack to compressed dense pack
  checker.command='gravityDeposition()'

def gravityDeposition():
 if O.iter==101:
  pred=pack.inCylinder((0,0,0),(0,0,H),Ro)-pack.inCylinder((0,0,0),(0,0,Hs),Rs)-pack.inCylinder((0,0,0),(0,0,H),Ri)
  assembly=pack.filterSpherePack(pred,sp,True)
  O.bodies.clear() #delete compressed pack sp

  print(len(O.bodies))
  assembly.toSimulation()
  print(len(O.bodies))
  for b in O.bodies:print(b.id,b.state.pos,b.shape.radius)
 else:
  print ('done')
  O.pause()
#...connect to next steps in simulation

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

> even though I could output their locations.

I could also output their position just after the toSimulation. However, if you print their positions the very next iteration before InsertionSortCollider, their positions are full of inf and nan values (I have no idea why, do some investigation).
Therefore you cannot see them (and also InsertionSortCollider "freezes").

cheers
Jan

Yuxuan Wang (ywang1) said : #6

Thanks you so much for the effort, Jan! I wasn't able to even locate the errors but you made it clear now. I'll keep researching for an answer on my end!

Best,
Yuxuan

Yuxuan Wang (ywang1) said : #7

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

Yuxuan Wang (ywang1) said : #8

Hi Jan,
Thank you again for the much needed help! I avoided the problem by outputting a table with the particles' locations and visualized the particles in a separate script.

Sincerely,
Yuxuan