How to reset the stiffness of the particles at a later stage?

Asked by 孙灿

Imagine that I am in a two-dimensional simulation box full of particles, and their materials are set, and the Young's modulus is also set, and we all know that the size of the Young's modulus affects the size of the stiffness, so the stiffness is also set at this time. Under dead weight conditions, these particles will pile together, after reaching equilibrium, I delete a circle range of particles between these particles, and then apply gravity again, this circular hole will collapse, but all I need is that this hole will deform without collapsing, can I change the particle stiffness (Young's modulus) around this hole so that this hole deforms without collapsing.

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
Jan Stránský (honzik) said :
#1

Hello,

> Imagine that I am in a two-dimensional simulation box full of particles

Why do you want us to imagine anything like that?

> can I change the particle stiffness (Young's modulus) around this hole

Yes.

Read [1] and provide relevant information (e.g. MWE) for more detailed answer (as it strongly depends on used materials, laws, etc.)

> so that this hole deforms without collapsing.

It will depend on many factors.
I am not sure if mere stiffness change can help here..

Cheers
Jan

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

Revision history for this message
孙灿 (suncan) said :
#2

#Material constants(1)
Density1 = 1836
FrictionAngle1 = 5
PoissonRatio1 = 0.4
Young1 = 1e8
Damp = 0.5
AvgRadius1 = 0.05
N_particles = 10000

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

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

SphereMat1 = O.materials.append(FrictMat(young = Young1, poisson = PoissonRatio1, frictionAngle = radians(FrictionAngle1), density = Density1))
WallMat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle)))

from yade import pack,plot
O.bodies.append(geom.facetBox((0.05, 15, 15), (0.1, 15, 15), wallMask=63,material=WallMat))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,30,30), rMean=0.1, rRelFuzz=0)
sp.toSimulation(material = SphereMat1)

(xdim,ydim,zdim)= aabbDim()
print("Height is ",zdim)

for b in O.bodies:
   if isinstance(b.shape,Sphere):
       b.state.blockedDOFs='ZxY'

       b.shape.color=(0,0,1.)

circleRadius=2
circleCenter = Vector3(0.05,15,6)
#myEngine.dead = True
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
        PyRunner(command='checkUnbalanced()', realPeriod=2,label="myEngine"),
        PyRunner(command='addPlotData()', iterPeriod=100),
 PyRunner(command='li()', realPeriod=2)
]
O.dt = 0.5 * PWaveTimeStep()
O.trackEnergy = True

def checkUnbalanced():
 if unbalancedForce() < .0001:
  O.pause()
                for b in O.bodies:
      if isinstance(b.shape,Sphere):
                        #b.state.blockedDOFs='zxy'
                        b.state.vel=(0,0,0)
                        b.state.angVel=(0,0,0)

def addPlotData():
 plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)

 for b in O.bodies:
  d = (b.state.pos - circleCenter).norm()
  if d < circleRadius:
   O.bodies.erase(b.id)
O.saveTmp()

from yade import qt
qt.Controller()
qt.View()

I just need this hole to have a strong stiffness so that it doesn't collapse immediately.

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

Please provide a MWE [1]:

W = working
IndentationError: unexpected indent
Please, try your code before posting. Especially that you posted code with the same problem several times before.

M = minimal, a few particles is sufficient to show how to change stiffness

Cheers
Jan

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

MWE atttached below.

Bodies:
======
to change body stiffness, assign new material.
(In general) do not do b.mat.young=newValue, since b.mat is (usually) shared among more particles and you would change young to them, too.

Interactions:
==========
Changing material has no effect on existing interactions.
One approach is to delete / erase existing interactions (completely or just those belonging to relevant particles). New interactions are created the very next step, reflecting new material.

This approach is not sufficient if you need "interaction history" (e.g. concerning friction).
Then you can use another approach - access interactions and change the values directly:
###
for i in body.intrs:
    i.phys.kn = newValue
###

MWE:
=====
####
mat1 = FrictMat(young=1)
mat2 = FrictMat(young=2)
spheres = s1,s2,s3,s4 = [sphere((x,0,0),.75,fixed=True,material=mat1) for x in (1,2,3,4)]

O.bodies.append(spheres)
O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    NewtonIntegrator(),
]
O.dt = 0

def printInfo():
    print("=============")
    print("step",O.iter)
    for i,body in enumerate(O.bodies):
        print("body",i,"young",body.mat.young)
    for i,intr in enumerate(O.interactions):
        print("interaction",intr.id1,intr.id2,"kn",intr.phys.kn)

O.step()
printInfo()

selectedSpheres = (s3,s4)
for body in selectedSpheres:
    body.mat = mat2
    for intr in body.intrs():
        O.interactions.erase(intr.id1,intr.id2)

O.step()
printInfo()
####

Cheers
Jan

Revision history for this message
孙灿 (suncan) said (last edit ):
#5

I know what you mean. But while I'm using this code, I don't know how to call the grain itself. The relevant code is as follows (not all of it):
shuzu=[]
circleRadius=2
circleCenter = Vector3(0.05,15,6)
for b in O.bodies:
 d = (b.state.pos - circleCenter).norm()
 if d < circleRadius:
  shuzu.append(?)
def printInfo():
    print("=============")
    print("step",O.iter)
    for i,body in enumerate(O.bodies):
        print("body",i,"young",body.mat.young)
    for i,intr in enumerate(O.interactions):
        print("interaction",intr.id1,intr.id2,"kn",intr.phys.kn)
O.step()
printInfo()
def li():
 for body in shuzu:
      body.mat = Mat2
      for intr in body.intrs():
          O.interactions.erase(intr.id1,intr.id2)
 printInfo()
I want to add eligible particles in the list shuzu[] using a for loop, and will change their material later using li(), I don't know how to call the particles themselves directly (put in the list shuzu[]), as I wrote shuzu.append(?). Maybe it's a very simple problem, but I don't know how to fix it.

In addition, can I achieve my purpose only through the operation you mentioned (make the hole have a certain degree of self-stability, so that the hole will slowly deform and finally collapse instead of collapsing quickly)

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

###
for b in O.bodies:
    ...
    shuzu.append(b)
###
?

Revision history for this message
孙灿 (suncan) said :
#7

shuzu.append(? ) means I'm not quite sure what to fill in () to add the particles themselves to the list. Or I add the id of the particle to the list, body.mat = Mat2 in the position of body I need to replace it with, otherwise an error will be reported in the line body.mat = Mat2, "'int' object has no attribute 'mat'"

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

> shuzu.append(? ) means I'm not quite sure what to fill in () to add the particles themselves to the list.

That is why I proposed to "add the particles themselves to the list" in #6.
Instant replay:
###
for b in O.bodies:
    ...
    shuzu.append(b)
###

> Or I add the id of the particle to the list

or, does not really matter.
Body and its ID is unique pair. Having ID you can easily get body itself with index access on body container: O.bodies[id]

> body.mat = Mat2 in the position of body I need to replace it with, otherwise an error will be reported in the line body.mat = Mat2, "'int' object has no attribute 'mat'"

sure, if body is ID (number), body.mat does not make sense

Cheers
Jan

Revision history for this message
孙灿 (suncan) said :
#9

shuzu=[]
circleRadius=2
circleCenter = Vector3(0.05,15,6)
for b in O.bodies:
 d = (b.state.pos - circleCenter).norm()
 if d < circleRadius:
  shuzu.append(b.id)
def printInfo():
    print("=============")
    print("step",O.iter)
    for i,body in enumerate(O.bodies):
        print("body",i,"young",body.mat.young)
    for i,intr in enumerate(O.interactions):
        print("interaction",intr.id1,intr.id2,"kn",intr.phys.kn)
O.step()
printInfo()
def li():
 for i in shuzu:
      O.bodies[i].mat = Mat2
      for intr in body.intrs():
          O.interactions.erase(intr.id1,intr.id2)
 printInfo()

I have changed the relevant code, but there is still the error 'Omega' object has no attribute 'bodise'.

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

> but there is still the error 'Omega' object has no attribute 'bodise'.

still? it is the first time you mention it.

It is simply expected error, as O.bodise does not exist (contrary to O.bodies)

Cheers
Jan

Revision history for this message
孙灿 (suncan) said :
#11

So how can I fix this?I'm not quite sure how to fix it.

Cheers

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

fix what?

'Omega' object has no attribute 'bodise'?
if so, simply do not use O.bodise

Cheers
Jan

Revision history for this message
孙灿 (suncan) said :
#13

>'Omega' object has no attribute 'bodise'?

What's the meaning of this? Does "Omega" here refer to the Greek letter Ω? What I add to the list is the particle ID (shuzu. append (b.id)), and I use O. bodies [I] Material=mat2, which cannot solve this problem.

The code is as follows:

#Material constants(1)
Density1 = 1836
FrictionAngle1 = 5
PoissonRatio1 = 0.4
Young1 = 1e8
Damp = 0.5
AvgRadius1 = 0.05
N_particles = 10000

#Material constants(23)
Density23 = 2000
FrictionAngle23 = 27
PoissonRatio23 = 0.35
Young23 = 1e9
Damp = 0.5
AvgRadius1 = 0.05
N_particles = 10000

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

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

mat1 = O.materials.append(FrictMat(young = Young1, poisson = PoissonRatio1, frictionAngle = radians(FrictionAngle1), density = Density1))
mat2 = O.materials.append(FrictMat(young = Young23, poisson = PoissonRatio1, frictionAngle = radians(FrictionAngle1), density = Density1))
WallMat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle)))

from yade import pack,plot
O.bodies.append(geom.facetBox((0.05, 15, 15), (0.5, 15, 15), wallMask=63,material=WallMat))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,30,30), rMean=0.3, rRelFuzz=0)
sp.toSimulation(material = mat1)

(xdim,ydim,zdim)= aabbDim()

print("Height is ",zdim)

for b in O.bodies:
   if isinstance(b.shape,Sphere):
       b.state.blockedDOFs='ZxY'
       b.shape.color=(0,0,1.)

circleRadius=2
circleCenter = Vector3(0.05,15,6)
circleRadius1=2.5
#myEngine.dead = True
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
        PyRunner(command='checkUnbalanced()', realPeriod=2,label="myEngine"),
        PyRunner(command='addPlotData()', iterPeriod=100)

]
O.dt = 0.5 * PWaveTimeStep()
O.trackEnergy = True

shuzu=[]
def checkUnbalanced():
    if unbalancedForce() < .0001:
        O.pause()
        zMax = max(b.state.pos[2] for b in O.bodies if isinstance(b.shape,Sphere))
        print("Z ",zMax)
        ceng = zMax-0.4
        for b in O.bodies:
            if isinstance(b.shape,Sphere):
                #b.state.blockedDOFs='zxy'
                b.state.vel=(0,0,0)
                b.state.angVel=(0,0,0)
            if b.state.pos[2]>ceng:
                print(b.id,b.state.pos)
                shuzu.append(b.id)
                #selectedSpheres+=(b,)
        print(shuzu)

   #zMax = max(b.state.pos[2] for b in O.bodies)

  #zMax = max(b.state.pos[2] for b in O.bodies)
  #print("Z ",zMax)
  #(xdim,ydim,zdim)= aabbDim()

  #print("Height is ",zdim)

def printInfo():
    print("=============")
    print("step",O.iter)
    for i,body in enumerate(O.bodies):
        print("body",i,"young",body.mat.young)
    for i,intr in enumerate(O.interactions):
        print("interaction",intr.id1,intr.id2,"kn",intr.phys.kn)
O.step()
printInfo()
def li():
    for i in shuzu:
        #i.mat = mat2
        O.bodies[i].material = mat2
    for intr in body.intrs():
        O.interactions.erase(intr.id1,intr.id2)
    printInfo()
def addPlotData():
    plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)

def duqu():
    for i in shuzu:
        b = O.bodies[i]
        print("voila",i,b.state.pos)
def wadong1():
    for b in O.bodies:
        d = (b.state.pos - circleCenter).norm()
        if d < circleRadius:
            O.bodies.erase(b.id)
O.saveTmp()

from yade import qt
qt.Controller()
qt.View()

The code is as follows: Error will be reported:
ArgumentError: Python argument types in
    BodyContainer.__getitem__(BodyContainer, Body)
did not match C++ signature:
    __getitem__(pyBodyContainer {lvalue}, int)

When I directly add the particles themselves to the list (shuzu. append (b)), and use i mat=mat2, error will still be reported:
ArgumentError: Python argument types in
    BodyContainer.__getitem__(BodyContainer, Body)
did not match C++ signature:
    __getitem__(pyBodyContainer {lvalue}, int)

The code is as follows:

#Material constants(1)
Density1 = 1836
FrictionAngle1 = 5
PoissonRatio1 = 0.4
Young1 = 1e8
Damp = 0.5
AvgRadius1 = 0.05
N_particles = 10000

#Material constants(23)
Density23 = 2000
FrictionAngle23 = 27
PoissonRatio23 = 0.35
Young23 = 1e9
Damp = 0.5
AvgRadius1 = 0.05
N_particles = 10000

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

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

mat1 = O.materials.append(FrictMat(young = Young1, poisson = PoissonRatio1, frictionAngle = radians(FrictionAngle1), density = Density1))
mat2 = O.materials.append(FrictMat(young = Young23, poisson = PoissonRatio1, frictionAngle = radians(FrictionAngle1), density = Density1))
WallMat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle)))

from yade import pack,plot
O.bodies.append(geom.facetBox((0.05, 15, 15), (0.5, 15, 15), wallMask=63,material=WallMat))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,30,30), rMean=0.3, rRelFuzz=0)
sp.toSimulation(material = mat1)

(xdim,ydim,zdim)= aabbDim()

print("Height is ",zdim)

for b in O.bodies:
   if isinstance(b.shape,Sphere):
       b.state.blockedDOFs='ZxY'
       b.shape.color=(0,0,1.)

circleRadius=2
circleCenter = Vector3(0.05,15,6)
circleRadius1=2.5
#myEngine.dead = True
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
        PyRunner(command='checkUnbalanced()', realPeriod=2,label="myEngine"),
        PyRunner(command='addPlotData()', iterPeriod=100)

]
O.dt = 0.5 * PWaveTimeStep()
O.trackEnergy = True

shuzu=[]
def checkUnbalanced():
    if unbalancedForce() < .0001:
        O.pause()
        zMax = max(b.state.pos[2] for b in O.bodies if isinstance(b.shape,Sphere))
        print("Z ",zMax)
        ceng = zMax-0.4
        for b in O.bodies:
            if isinstance(b.shape,Sphere):
                #b.state.blockedDOFs='zxy'
                b.state.vel=(0,0,0)
                b.state.angVel=(0,0,0)
            if b.state.pos[2]>ceng:
                print(b.id,b.state.pos)
                shuzu.append(b)
                #selectedSpheres+=(b,)
        print(shuzu)

   #zMax = max(b.state.pos[2] for b in O.bodies)

  #zMax = max(b.state.pos[2] for b in O.bodies)
  #print("Z ",zMax)
  #(xdim,ydim,zdim)= aabbDim()

  #print("Height is ",zdim)

def printInfo():
    print("=============")
    print("step",O.iter)
    for i,body in enumerate(O.bodies):
        print("body",i,"young",body.mat.young)
    for i,intr in enumerate(O.interactions):
        print("interaction",intr.id1,intr.id2,"kn",intr.phys.kn)
O.step()
printInfo()
def li():
    for i in shuzu:
        i.mat = mat2
        #O.bodies[i].material = mat2
    for intr in body.intrs():
        O.interactions.erase(intr.id1,intr.id2)
    printInfo()
def addPlotData():
    plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)

def duqu():
    for i in shuzu:
        b = O.bodies[i]
        print("voila",i,b.state.pos)
def wadong1():
    for b in O.bodies:
        d = (b.state.pos - circleCenter).norm()
        if d < circleRadius:
            O.bodies.erase(b.id)
O.saveTmp()

from yade import qt
qt.Controller()
qt.View()

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

>> 'Omega' object has no attribute 'bodise'?
> What's the meaning of this?

exactly as the error says, 'Omega' object has no attribute 'bodise'.
You try to access "bodise" attribute of Omega object, but Omega object has no attribute "bodise"

> Does "Omega" here refer to the Greek letter Ω?

No.
Omega is Yade class [2], basically "The simulation". Usually Omega is used as "O" object in Yade scripts, (O.bodies, O.interactions, O.engines, ... ).
Try
print(O)
in Yade terminal.

> What I add to the list is the particle ID (shuzu. append (b.id)), and I use O. bodies [I] Material=mat2, which cannot solve this problem.

This is invalid syntax, you cannot expect it to solve any problem..

>> 'Omega' object has no attribute 'bodise'?
> The code is as follows:

you have no "bodise" in the code, please be consistent

> Error will be reported:
> ArgumentError: Python argument types in
> BodyContainer.__getitem__(BodyContainer, Body)
> did not match C++ signature:
> __getitem__(pyBodyContainer {lvalue}, int)

It says, that O.bodies[arg] expects "arg" to be int, but actually it is Body.
To solve it, put arg.id (Body.id as int)

Why do you mention the same error twice?

> and use i mat=mat2, error will still be reported: ArgumentError

if you use this, SyntaxError should be reported.

Moreover, I have no problem with both codes you have posted..

Cheers
Jan

[2] https://yade-dem.org/doc/yade.wrapper.html#omega

Revision history for this message
孙灿 (suncan) said :
#15

>It says, that O.bodies[arg] expects "arg" to be int, but actually it is Body.
>To solve it, put arg.id (Body.id as int)

I try to output b The type of id is found to be 'int' instead of 'body'. When using the li() function, I output the type of the body in for body in shuzu: again, which is still 'int'. Why does' that O. bodies [arg] expect "arg" to be int, but actually it is Body. 'still occur.

>Moreover, I have no problem with both codes you have posted..

In the previous code, I put the functions outside checkUnbalanced() (I manually input the functions). Because the error is in li(), you can't find an error by clicking Run. Now I have readjusted the code as follows:

#Material constants(1)
Density1 = 1836
FrictionAngle1 = 5
PoissonRatio1 = 0.4
Young1 = 1e8
Damp = 0.5
AvgRadius1 = 0.05
N_particles = 10000

#Material constants(23)
Density23 = 2000
FrictionAngle23 = 27
PoissonRatio23 = 0.35
Young23 = 1e9
Damp = 0.5
AvgRadius1 = 0.05
N_particles = 10000

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

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

mat1 = O.materials.append(FrictMat(young = Young1, poisson = PoissonRatio1, frictionAngle = radians(FrictionAngle1), density = Density1))
mat2 = O.materials.append(FrictMat(young = Young23, poisson = PoissonRatio1, frictionAngle = radians(FrictionAngle1), density = Density1))
WallMat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle)))

from yade import pack,plot
O.bodies.append(geom.facetBox((0.05, 15, 15), (0.5, 15, 15), wallMask=63,material=WallMat))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,30,30), rMean=0.3, rRelFuzz=0)
sp.toSimulation(material = mat1)

(xdim,ydim,zdim)= aabbDim()

print("Height is ",zdim)

for b in O.bodies:
   if isinstance(b.shape,Sphere):
       b.state.blockedDOFs='ZxY'
       b.shape.color=(0,0,1.)

circleRadius=2
circleCenter = Vector3(0.05,15,6)
circleRadius1=2.5
#myEngine.dead = True
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
        PyRunner(command='checkUnbalanced()', realPeriod=2,label="myEngine"),
        PyRunner(command='addPlotData()', iterPeriod=100)

]
O.dt = 0.5 * PWaveTimeStep()
O.trackEnergy = True

shuzu=[]
def checkUnbalanced():
    if unbalancedForce() < .0001:
        O.pause()
        zMax = max(b.state.pos[2] for b in O.bodies if isinstance(b.shape,Sphere))
        print("Z ",zMax)
        ceng = zMax-0.4
        for b in O.bodies:
            if isinstance(b.shape,Sphere):
                #b.state.blockedDOFs='zxy'
                b.state.vel=(0,0,0)
                b.state.angVel=(0,0,0)
            if b.state.pos[2]>ceng:
                print(b.id,b.state.pos)
                #a=int(b.id)
                #print(type(b.id))
                shuzu.append(b.id)
  print(type(b.id))
                #selectedSpheres+=(b,)
        print(shuzu)

   #zMax = max(b.state.pos[2] for b in O.bodies)

  #zMax = max(b.state.pos[2] for b in O.bodies)
  #print("Z ",zMax)
  #(xdim,ydim,zdim)= aabbDim()

  #print("Height is ",zdim)
        li()

def printInfo():
    print("=============")
    print("step",O.iter)
    for i,body in enumerate(O.bodies):
        print("body",i,"young",body.mat.young)
    for i,intr in enumerate(O.interactions):
        print("interaction",intr.id1,intr.id2,"kn",intr.phys.kn)
O.step()
printInfo()
def li():
    for body in shuzu:
        print(type(body))
        #body.mat = mat2
        O.bodies[body].material = mat2
    for intr in body.intrs():
        O.interactions.erase(intr.id1,intr.id2)
    printInfo()
def addPlotData():
    plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)

def duqu():
    for i in shuzu:
        b = O.bodies[i]
        print("voila",i,b.state.pos)
def wadong1():
    for b in O.bodies:
        d = (b.state.pos - circleCenter).norm()
        if d < circleRadius:
            O.bodies.erase(b.id)
O.saveTmp()

from yade import qt
qt.Controller()
qt.View()

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

> Now I have readjusted the code as follows:

IndentationError: unindent does not match any outer indentation level

> Why does' that O. bodies [arg] expect "arg" to be int, but actually it is Body. 'still occur.

Usage of O.bodies[...] seems perfectly OK in the provided code.
Isn't the error different? like "AttributeError: 'int' object has no attribute 'intrs'"?

> you can't find an error by clicking Run

please always try to provide a MWE, W = working, reproducing your problem..

> def li():
> for body in shuzu:
> print(type(body))
> #body.mat = mat2
> O.bodies[body].material = mat2
> for intr in body.intrs():
> O.interactions.erase(intr.id1,intr.id2)

compare your current code (levels of the for loops) with #4

Cheers
Jan

Revision history for this message
孙灿 (suncan) said (last edit ):
#17

Sorry, I know why my published code always has indentation errors.
The code is as follows (it should work this time):

#Material constants(1)
Density1 = 1836
FrictionAngle1 = 5
PoissonRatio1 = 0.4
Young1 = 1e8
Damp = 0.5
AvgRadius1 = 0.05
N_particles = 10000

#Material constants(23)
Density23 = 2000
FrictionAngle23 = 27
PoissonRatio23 = 0.35
Young23 = 1e9
Damp = 0.5
AvgRadius1 = 0.05
N_particles = 10000

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

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

mat1 = O.materials.append(FrictMat(young = Young1, poisson = PoissonRatio1, frictionAngle = radians(FrictionAngle1), density = Density1))
mat2 = O.materials.append(FrictMat(young = Young23, poisson = PoissonRatio1, frictionAngle = radians(FrictionAngle1), density = Density1))
WallMat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle)))

from yade import pack,plot
O.bodies.append(geom.facetBox((0.05, 15, 15), (0.5, 15, 15), wallMask=63,material=WallMat))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0.05,0,0),Vector3(0.05,30,30), rMean=0.3, rRelFuzz=0)
sp.toSimulation(material = mat1)

(xdim,ydim,zdim)= aabbDim()

print("Height is ",zdim)

for b in O.bodies:
   if isinstance(b.shape,Sphere):
       b.state.blockedDOFs='ZxY'
       b.shape.color=(0,0,1.)

circleRadius=2
circleCenter = Vector3(0.05,15,6)
circleRadius1=2.5
#myEngine.dead = True
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
        PyRunner(command='checkUnbalanced()', realPeriod=2,label="myEngine"),
        PyRunner(command='addPlotData()', iterPeriod=100)

]
O.dt = 0.5 * PWaveTimeStep()
O.trackEnergy = True

shuzu=[]
def checkUnbalanced():
    if unbalancedForce() < .0001:
        O.pause()
        zMax = max(b.state.pos[2] for b in O.bodies if isinstance(b.shape,Sphere))
        print("Z ",zMax)
        ceng = zMax-0.4
        for b in O.bodies:
            if isinstance(b.shape,Sphere):
                #b.state.blockedDOFs='zxy'
                b.state.vel=(0,0,0)
                b.state.angVel=(0,0,0)
            if b.state.pos[2]>ceng:
                print(b.id,b.state.pos)
                #a=int(b.id)
                #print(type(b.id))
                shuzu.append(b.id)
                print(type(b.id))
                #selectedSpheres+=(b,)
        print(shuzu)
        li()

def printInfo():
    print("=============")
    print("step",O.iter)
    for i,body in enumerate(O.bodies):
        print("body",i,"young",body.mat.young)
    for i,intr in enumerate(O.interactions):
        print("interaction",intr.id1,intr.id2,"kn",intr.phys.kn)
O.step()
printInfo()
def li():
    for body in shuzu:
        print(type(body))
        #body.mat = mat2
        O.bodies[body].material = mat2
        for intr in O.bodies[body].intrs():
            O.interactions.erase(intr.id1,intr.id2)
    printInfo()
def addPlotData():
    plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)
O.saveTmp()

from yade import qt
qt.Controller()
qt.View()

>Isn't the error different? like "AttributeError: 'int' object has no attribute 'intrs'"?
No, it will still report errors mentioned above.

>compare your current code (levels of the for loops) with #4
Yes, I found an error (for intr in body. intrs():) with indentation error. I changed it in the code mentioned above.

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

> No, it will still report errors mentioned above.

No, it will report a different error:

ArgumentError: Python argument types in
    None.None(Body, int)
did not match C++ signature:
    None(yade::Body {lvalue}, boost::shared_ptr<yade::Material>)

Cheers
Jan

Revision history for this message
孙灿 (suncan) said :
#19

>ArgumentError: Python argument types in
 >None.None(Body, int)
>did not match C++ signature:
 >None(yade::Body {lvalue}, boost::shared_ptr<yade::Material>)

Yes, that's the error, I don't have a way to solve it.

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

> No, it will still report errors mentioned above.
> Yes, that's the error

It is a good habit to read error messages.
Even if you do not fully understand them, you should be able to distinguish that it is different error.

> I don't have a way to solve it.

The error gives you exact line where the problem is.
Even if you do not fully understand the error, there are some clues making it possible for you to solve it yourself with just a little debugging.

The keywords of the error message are "Body", "int", "Material".
Combining this info with the error line
O.bodies[body].material = mat2
it is not difficult after some debugging (a few print(...) is sufficient) to realize that mat2 is int.

In an independent Yade session you can test something like:
###
b = sphere((0,0,0),1)
print(b.material)
b.material = FrictMat() # OK
b.material = 1 # the same error message
###

Similarly to bodies, material and its ID is unique pair.
Similarly to O.bodies.append, O.material.append returns material ID, not material itself

So either use this pattern for materials
###
mat = SomeMat(...)
matID = O.materials.append(mat)
###
or
###
matID = O.materials.append(...)
mat = O.materials[matID]
###

Cheers
Jan

Revision history for this message
孙灿 (suncan) said :
#21

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