Can not see existed particles in "show 3D"

Asked by Leonard on 2019-09-10

Hi,
I'm trying to generate a cylinder sample with target porosity, the MWE comes from [1]. The main changes are:
1. changed the sphere material from FrictMat to CohFrictMat, as well as corresponding engines,
2. cut the cubic sample to fulfil a cylinder shape after target porosity is reached.
The problem is after the simulation finished (run about 30 seconds), I can not find any particles in the scene, while I use 'len(O.bodies)' to check, I can get that there are 870 bodies. Here is the MWE:
########
from yade import pack
num_spheres=1000
key='_triax_base_'
targetPorosity = 0.43
compFricDegree = 30
finalFricDegree = 30
rate=-0.02
damp=0.2
stabilityThreshold=0.01
young=5e6
mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0.14)

cement = CohFrictMat(young=30e9,poisson=0.3,frictionAngle=radians(30),density=2650.0,isCohesive=True,normalCohesion=1e8, shearCohesion=1e8,label='cement')
walls = FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls')

O.materials.append(cement)
O.materials.append(walls)

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()

x=[5.0,6.0005088055510125,6.807147320174279,7.200507128170075,7.399996645238125,7.604595185078019,7.809172757656196,8.008676252898724,8.20306373628217,8.41275033163219,9.005096442414814]
y=[0,9.972588799847912,20.08747541588565,29.9397121334748,39.93424666725376,50.065627529175956,59.78605006143409,70.0545568149891,80.04923113051282,89.9064999909142,100.0]

#### change the size units
psdSizes=x
for i in range(0,len(x)):
    psdSizes[i]=x[i]*0.001

## change the unit to percentage
psdCumm=y
for i in range(0,len(y)):
    psdCumm[i]=y[i]*0.01

sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,seed=1)
sp.toSimulation(material='cement')

Gl1_Sphere.quality=3

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=False,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(),Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key),
 newton
]
triax.goal1=triax.goal2=triax.goal3=-10000

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
    break

print "### Isotropic state saved ###"

import sys
while triax.porosity>targetPorosity:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r compFrictionDegree: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)

print "### Compacted state saved ###"

aabb=utils.aabbExtrema()
predX=(aabb[0][0]+aabb[1][0])/2
predY=(aabb[0][1]+aabb[1][1])/2
predR=min((aabb[1][0]-aabb[0][0])/2,(aabb[1][1]-aabb[0][1])/2)

pred=pack.inCylinder((predX,predY,0),(predX,predY,0.14),predR)

## erase Walls and spheres which is not in cylinder
def cylinderShape():
 for i in O.bodies:
  if not isinstance(i.shape,Sphere):
   O.bodies.erase(i.id)
  if isinstance(i.shape,Sphere):
    if not pred(i.state.pos, i.shape.radius): # if sphere with c,r is in pred(cylinder), then add this sphere into bodies.
     O.bodies.erase(i.id)
cylinderShape()
########

Could you please help me with this problem?
Thanks,
Leonard

[1]https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Leonard
Solved:
2019-09-11
Last query:
2019-09-11
Last reply:
2019-09-10
Jan Stránský (honzik) said : #1

what is the AABB of the final state compared to the initial state? [1,2]
Jan

[1] https://yade-dem.org/doc/yade.utils.html#yade.utils.aabbDim
[2] https://yade-dem.org/doc/yade.utils.html#yade._utils.aabbExtrema

Leonard (z2521899293) said : #2

Hi Jan,
Thanks for your reply.
In my case, I just use AABB to provide parameters for creating a cylinder shape, you may ignore it, let's comment the last part, as the MWE below. However. I meet the same problem (Can not see existed particles in "show 3D" ).
########
from yade import pack
num_spheres=1000
key='_triax_base_'
targetPorosity = 0.43
compFricDegree = 30
finalFricDegree = 30
rate=-0.02
damp=0.2
stabilityThreshold=0.01
young=5e6
mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0.14)

cement = CohFrictMat(young=30e9,poisson=0.3,frictionAngle=radians(30),density=2650.0,isCohesive=True,normalCohesion=1e8, shearCohesion=1e8,label='cement')
walls = FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls')

O.materials.append(cement)
O.materials.append(walls)

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()

x=[5.0,6.0005088055510125,6.807147320174279,7.200507128170075,7.399996645238125,7.604595185078019,7.809172757656196,8.008676252898724,8.20306373628217,8.41275033163219,9.005096442414814]
y=[0,9.972588799847912,20.08747541588565,29.9397121334748,39.93424666725376,50.065627529175956,59.78605006143409,70.0545568149891,80.04923113051282,89.9064999909142,100.0]

#### change the size units
psdSizes=x
for i in range(0,len(x)):
    psdSizes[i]=x[i]*0.001

## change the unit to percentage
psdCumm=y
for i in range(0,len(y)):
    psdCumm[i]=y[i]*0.01

sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,seed=1)
sp.toSimulation(material='cement')

Gl1_Sphere.quality=3

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=False,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(),Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key),
 newton
]
triax.goal1=triax.goal2=triax.goal3=-10000

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
    break

print "### Isotropic state saved ###"

import sys
while triax.porosity>targetPorosity:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r compFrictionDegree: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)

print "### Compacted state saved ###"
############

Leonard

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

Hi Leonard,
sorry for not very clear explanation. What I meant is that particles after simulation may be "expolded", far away from each other, in 3D view seemingly not there (too small, the view tries to fit all scene). Aabb values after the simulation can (dis)prove this reason.

> The problem is after the simulation finished (run about 30 seconds), I can not find any particles

do you see particles before the simulation starts?

Jan

Leonard (z2521899293) said : #4

Hi Jan,
Sorry for my unclear explanation before.
The initial AABB is:
(Vector3(0.00005412909461737815,7.493929978427943e-6,0.000025167230754637605),
 Vector3(0.0699463745256682,0.06994859708381397,0.13982812911732453))
The final AABB is:
(Vector3(-5975.64011320338,-8044.161847180958,-7929.333516594808),
 Vector3(6450.074621608338,8773.927665611896,6573.516040778314))
I think you are right that the particles are exploded.
Could you give some advice how to solve this problem?
Thanks
Leonard

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

> Could you give some advice how to solve this problem?

not really, sorry :-)
open a new question mentioning the actual problem of "explosion".

cheers
Jan

Leonard (z2521899293) said : #6

Thanks Jan,
I reviewed some good answers and find two ways to solve it:
1. increase wall's stiffness ;
2. reduce O.dt.

For this MWE, increasing stiffness could directly work.
Thank you again for helping me diagnose the script.

cheers
Leonard

Leonard (z2521899293) said : #7

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