randomdensepack exclusion from walls

Asked by JOHN on 2017-12-06

Good evening,
Previously, i had asked how to create a group of particles and exclude those that are inside the walls.
This was the solution

  pred = inGtsSurface(s)
  sp=pack.SpherePack()
  sp.makeCloud((10,40,3),(30,60,12),rMean=0.5)

  # remove spheres completely inside walls
  for c,r in sp:
     if pred(c,0):
        continue
     O.bodies.append(sphere(c,r))

  # remove spheres partially inside walls
  O.dt = 0
  O.step() # interactions are created afterwards
  toErase = set()
  for i in O.interactions:
     b1,b2 = [O.bodies[i] for i in (i.id1,i.id2)]
     if any(isinstance(b.shape,Facet) for b in (b1,b2)): # if facet is involved, delete
        toErase.add(b1)
        toErase.add(b2)
  toErase = [b for b in toErase if isinstance(b.shape,Sphere)] # delete just spheres
  for b in toErase: # delete the spheres
  O.bodies.erase(b.id)
  O.engines=[

However, a much denser particle cloud is needed
To that effect i subsituted

  sp=pack.SpherePack()
  sp.makeCloud((10,40,3),(30,60,12),rMean=0.5)

with

  sp=pack.randomDensePack(pack.inAlignedBox((-10,40,3),(20,50,10)),spheresInCell=1000,radius=.05,memoizeDb='/tmp/triaxPackCache.sqlite')

ln: failed to create symbolic link 'yadeimport.py': File exists
False
True
Found suitable packing in /tmp/triaxPackCache.sqlite (radius=0.05±0,N=1000,dim=2.2861×0.762034×0.533423,periodic,scale=1), created Wed Dec 6 15:34:07 2017
/usr/lib/x86_64-linux-gnu/yade/py/yade/pack.py:296: FutureWarning: The default behavior will change; specify returnSpherePack=True for the new behavior, and False to get rid of this warning (your code will break in the future, however). The returned SpherePack object can be added to the simulation using SpherePack.toSimulation()
  warnings.warn('The default behavior will change; specify returnSpherePack=True for the new behavior, and False to get rid of this warning (your code will break in the future, however). The returned SpherePack object can be added to the simulation using SpherePack.toSimulation()',category=FutureWarning)
Killed

Any help is appreciated
Thank you for your time
John

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
JOHN
Solved:
2017-12-06
Last query:
2017-12-06
Last reply:
2017-12-06
Jan Stránský (honzik) said : #1

Hi John,
please specify what is actually the problem
- ln: failed ... ?
- warnings of randomDensePack?
- Killed ?

just a note, in the first case you use radius 0.5, while in the second case 0.05, which is roughly (including makeCloud x randomDensePack effect) 5000 times more particles..

cheers
Jan

JOHN (washingmachine) said : #2

Thank you for the reply,
the problem is that it doesnt output anything
i changed the radius to
    sp=pack.randomDensePack(pack.inAlignedBox((-10,40,3),20,50,10)),spheresInCell=1000,radius=.5,memoizeDb='/tmp/triaxPackCache.sqlite')

and now i get the much more informative

Traceback (most recent call last):
  File "calltest.py", line 8, in <module>
    run()
  File "/home/john/Desktop/yade/DONE/test.py", line 44, in run
    for c,r in sp:
TypeError: 'Body' object is not iterable

(+the rest of the warnigns)

my full code is

import sys
import os

os.system("ln -s /usr/bin/yade yadeimport.py")
from yadeimport import yade

import gts
from yadeimport import *
from yade import ymport
from yade import qt
import math
def printing():
 for b in O.bodies:
  if isinstance(b.shape,Sphere):
        print b.id,b.state.pos[0],b.state.pos[1],b.state.pos[2]

def run():

  facets = ymport.stl('maze4lid.stl')
  rod1 = O.bodies.append(facets)

  # converts facets to gts (see the other question)
  s = gts.Surface()
  for facet in facets:
     vs = [facet.state.pos + facet.state.ori*v for v in facet.shape.vertices]
     vs = [gts.Vertex(v[0],v[1],v[2]) for v in vs]
     es = [gts.Edge(vs[i],vs[j]) for i,j in ((0,1),(1,2),(2,0))]
     f = gts.Face(es[0],es[1],es[2])
     s.add(f)
  print s.is_closed()
  threshold = 1e-3
  s.cleanup(threshold)
  print s.is_closed()
  assert s.is_closed()
  a=(0,-9.81,0)
  # use gts to filter spheres
  pred = inGtsSurface(s)
    sp=pack.randomDensePack(pack.inAlignedBox((-10,40,3),(20,50,10)),spheresInCell=1000,radius=.5,memoizeDb='/tmp/triaxPackCache.sqlite')
  # remove spheres completely inside walls
  for c,r in sp:
     if pred(c,0):
        continue
     O.bodies.append(sphere(c,r))

  # remove spheres partially inside walls
  O.dt = 0
  O.step() # interactions are created afterwards
  toErase = set()
  for i in O.interactions:
     b1,b2 = [O.bodies[i] for i in (i.id1,i.id2)]
     if any(isinstance(b.shape,Facet) for b in (b1,b2)): # if facet is involved, delete
        toErase.add(b1)
        toErase.add(b2)
  toErase = [b for b in toErase if isinstance(b.shape,Sphere)] # delete just spheres
  for b in toErase: # delete the spheres
     O.bodies.erase(b.id)
  O.engines=[

     ForceResetter(),
     InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
     InteractionLoop(
  # handle sphere+sphere and facet+sphere collisions
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
     ),
     NewtonIntegrator(gravity=(0,-9.81,0),damping=0.4, label='newtonInt'),
    # PyRunner(iterPeriod=5,command='printing()'),
      VTKRecorder(fileName='/home/john/Desktop/yade/DONE/vis/3d-vtk-',recorders=['spheres'],iterPeriod=100),
        PyRunner(iterPeriod=1,command='gravity(a)'),
      PyRunner(iterPeriod=1,command='checkAndDelete()'),

  ]
  O.dt=.8*PWaveTimeStep()
  O.step()

def irun(num):

  O.run(num,1)
  print O.time
  #printing()

def strategy():
  a=(9.81,0,0)
  return a
def gravector(a):
 b=a

def gravity(a):

   newtonInt.gravity = a

def checkAndDelete():
    for b in O.bodies:
  if isinstance(b.shape,Sphere):
         x,y,z = b.state.pos
        if y<-90: # modify this condition to your needs
            O.bodies.erase(b.id)

and is called via a python script

import sys
import os
from test import *
from yadeimport import *

a=(9.81,0,0)
run()
irun(1000)
a=(0,-9.81,0)
irun(3000)

(which worked untill the randomdensepack addition)

Cheers
John

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

use (as one of the warning says)
sp=pack.randomDensePack(..., returnSpherePack=True)
cheers
Jan

JOHN (washingmachine) said : #4

Thank you!
apparently it also died because the 0.05 radius was too much for my machine to handle so
double thank you
cheers
john