How to generate a clump

Asked by William

Hi,

The program doesn't work but I can't find out what's wrong.

#Wall constants
WDensity = 0
WFrictionAngle = 0.0
WPoissonRatio = 0.5
WYoung = 50e9
Mat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle), density = WDensity))

axis = (0,0,1)
angle = 0.25*pi
topbodylist = []
for i in range(10):
    box = utils.box((1+2*i,1,1),(sqrt(2)/2,sqrt(2)/2,sqrt(2)/2))
    box.mat=O.materials[Mat]
    box.state.ori = axis,angle
    boxid = O.bodies.append(box)
    topbodylist.append(boxid)
topwallclump = O.bodies.clump(topbodylist)

Question information

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

just run the following code, it works.

axis = (0,0,1)
angle = 0.25*pi
topbodylist = []
for i in range(10):
    box = utils.box((1+2*i,1,1),(sqrt(2)/2,sqrt(2)/2,sqrt(2)/2))
    box.state.ori = axis,angle
    boxid = O.bodies.append(box)
    topbodylist.append(boxid)
topwallclump = O.bodies.clump(topbodylist)

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

Hello,

if you get an error, please provide that error ([1], point 2).

I have tried your script with

python3.6: /builds/yade-dev/trunk/deb/yadedaily/core/Clump.cpp:241: static void yade::Clump::updateProperties(const boost::shared_ptr<yade::Body>&, unsigned int): Assertion `M > 0' failed.

even being a bit cryptic, "Assertion `M > 0' failed" means that "M" is either 0 or negative, but should be positive (and non-zero).

> just run the following code, it works.

what is different from the erroneous code? Material.

Together with "M > 0 failed" and
WDensity = 0
one can deduce, that "M" is mass (zero for zero density).

Either set mass of boxes manually non-zero or use non-zero value for density and your code works.

Cheers
Jan

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

Revision history for this message
William (qfxx-123) said :
#3

Hi, Jan
My error is same with yours.

>Either set mass of boxes manually non-zero or use non-zero value for density and your code works.

It doesn't work.

I have two scripts and they should run in sequence. The first one can be run successfully, but the 2nd one still doesn't work due to above problem.

the first script:

#Material constants
Density = 3000
FrictionAngle = 45
PoissonRatio = 0.5
Young = 300e6
Damp = 0.5
AvgRadius = 0.00004
N_particles = 12000

#Wall constants
WDensity = 0
WFrictionAngle = 0.0
WPoissonRatio = 0.5
WYoung = 50e9

#Packing variables
mn = Vector3(0,0.001,0)
mx = Vector3(0.1,0.041,0.006)

#Confining variables
ConfPress1 = -0.2e5 #pre-compression
ConfPress = -1.0e5

#time calculation
startT = O.time
endT = O.time
timeSpent = endT - startT

########################################
#import necessary packages
from yade import pack,plot,os,timing
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab

########################################
### Sample Preparing ###################
########################################

#Create materials for spheres and plates
SphereMat = O.materials.append(FrictMat(young = Young, poisson = PoissonRatio, frictionAngle = radians(FrictionAngle), density = Density))
WallMat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle), density = WDensity))

#Create walls for packing
wallIds = O.bodies.append(aabbWalls([mn,mx],thickness=0.000001,material=WallMat))

sp = pack.SpherePack()
psdSizes,psdCumm=[0.000059,0.00006,0.00023,0.00024],[0,0.5,0.5,1]
sp.makeCloud(Vector3(0,0.001,0.003),Vector3(0.1,0.041,0.003),num=N_particles,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True)
sp.toSimulation(material = SphereMat)

O.usesTimeStepper=True
O.trackEnerty=True
#################################
#####Defining triaxil engines####
#################################

###first step: compression#######
triax1=TriaxialStressController(
    wall_back_activated = False, #for 2d simulation
    wall_front_activated = False,
    thickness = 0.000001,
    maxMultiplier=1.+1.5e5/Young, # spheres growing factor (fast growth)
    finalMaxMultiplier=1.+4e3/Young,
    internalCompaction = True, # If true the confining pressure is generated by growing particles
    stressMask = 7,
    computeStressStrainInterval = 10,

    goal1 = ConfPress1,
    goal2 = ConfPress1,
)

newton=NewtonIntegrator(damping=Damp)

###engine
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()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
    triax1,
    newton,
    PyRunner(realPeriod=10,command='checkUnbalanced()',label='check'),
]
# Simulation stop conditions defination
def checkUnbalanced():
    unb=unbalancedForce()
    mStress = (triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2.0
    s1 = triax1.stress(triax1.wall_right_id)[0]
    s2 = triax1.stress(triax1.wall_top_id)[1]
    print(unb,s2,s1)
    if unb<0.01 and abs(ConfPress1-s2)/abs(ConfPress1)<0.05:
       O.pause()
       O.save('initial.yade.bz2')

the second script:

#Wall constants
WDensity = 3000
WFrictionAngle = 0.001
WPoissonRatio = 0.5
WYoung = 50e9
Mat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle), density = WDensity))

#restart
O.load('initial.yade.bz2')

for b in range(6):
    O.bodies.erase(b)

axis = (0,0,1)
angle = 0.25*pi
bottombodylist = []
for i in range(100):
    box = utils.box((0.0005+0.001*i,.0005,0.003),(0.0005*sqrt(2)/2,0.0005*sqrt(2)/2,0.0005*sqrt(2)/2))
    box.mat=O.materials[Mat]
    box.state.ori = axis,angle
    boxid = O.bodies.append(box)
    bottombodylist.append(boxid)
bottomwallclump = O.bodies.clump(bottombodylist)

for b in O.bodies:
    if isinstance(b.shape,Sphere):
        b.shape.color=(0.,0.,1.)
for i in range(-100,-1):
    O.bodies[i].shape.color=(1.,0.,1.)

Revision history for this message
William (qfxx-123) said :
#4

If delete "bottomwallclump = O.bodies.clump(bottombodylist)", it works, but I need them to be a clump.

In an other way to ask, how could I turn these boxes into a clump?

the second script which works:

#Wall constants
WDensity = 3000
WFrictionAngle = 0.001
WPoissonRatio = 0.5
WYoung = 50e9
Mat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle), density = WDensity))

#restart
O.load('initial.yade.bz2')

for b in range(6):
    O.bodies.erase(b)

axis = (0,0,1)
angle = 0.25*pi
bottombodylist = []
for i in range(100):
    box = utils.box((0.0005+0.001*i,.0005,0.003),(0.0005*sqrt(2)/2,0.0005*sqrt(2)/2,0.0005*sqrt(2)/2))
    box.mat=O.materials[Mat]
    box.state.ori = axis,angle
    boxid = O.bodies.append(box)
    bottombodylist.append(boxid)

for b in O.bodies:
    if isinstance(b.shape,Sphere):
        b.shape.color=(0.,0.,1.)
for i in range(-100,-1):
    O.bodies[i].shape.color=(1.,0.,1.)

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

> It doesn't work.

you have to use is correctly, i.e. actually using the material with non-zero density

> box.mat=O.materials[Mat]

this has no effect on already set values, e.g. box.state.mass

> O.load('initial.yade.bz2')

previously appended
Mat = O.materials.append(FrictMat(...))
is overwritten with loaded materials, i.e.
SphereMat = O.materials.append(FrictMat(..., density = Density))
WallMat = O.materials.append(FrictMat(..., density = WDensity)) # 0 !!

> box = utils.box(...)

if material is not explicitly specified, lastly appended material is used.
In this case, WallMat with zero density
Therefore the boxes have zero mass and you get the error.

Solution: e.g. append Mat after O.load (such that it is the lastly appended material)

Cheers
Jan

Revision history for this message
William (qfxx-123) said :
#6

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