intRadius(aabbenlargefactor) doesn't work

Hi,

Unfortunately, your question / script is way too long to make it easy for people (me, at least) to help easily. It's good to provide some details, but there are too many here, in my opinion (you may give a look to https://www.yade-dem.org/wiki/Howtoask to have some ideas how to improve)

Anyway, I think aabbenlargefactor "works" but there is most probably a problem in your script (I could not check because of length). As a first guess, I would like to remind you you have to run some DEM iterations (e.g. O.step()) while "intRadius" > 1 for interactions to be detected with this interaction radius.

Also, I would recommend to try first with a much simpler/shorter script with just two non-touching spheres, and play with aabbEnlargeFactor to nevertheless create the interaction.

Thanks Jérôme Duriez for you kind advice.
I will edit the question.

Hello,

thanks for reducing the script, but next time please do it in more
sophisticated way. Your question is about intRadius, but in your script
there is nothing like this.. :-)

There are a few problems in your scripts:

As Jerome pointed, there should be one step between setting intRadius and
counting interactions (new interactions need collider to run and
InteractionLoop to decide which interactions are real and which only
potential etc.). Just instant setting of a parameter does not creates new
interactions by itself. So number of interactions in isotropic_compaction()
function is not "before reset intR", but actually before setting it :-)

So you will get actual number of interactions in set_interactions()
function, but the number will be before reset of intR (as there is not step
between interaction counting and intR resetting) :-)

I tried your original script with the behavior you described. However, if I
set normalCohesion to 4500e10, I got approx double interactions. I am not
familiar with CohFrictMat, so you will have to ask somebody else how to
deal with this problem.

As Jerome proposed, a small script (not only here, but in general) could be
a way to understand how things works..

cheers
Jan

2016-11-17 4:43 GMT+01:00 weijy <email address hidden>:

> Question #404135 on Yade changed:
> https://answers.launchpad.net/yade/+question/404135
>
> Description changed to:
> Hi All,
>
> I am a beginner and have been trying to remodel the triaxial compression
> process in below paper:
> [1] Scholt E S L and Donz E F E D E. 2013. A DEM model for soft and hard
> rocks: Role of grain interlocking on strength [J]. J. Mech. Phys.
> Solids.,61(2): 352-369.
> In chapter 2.2 interaction range is introduced.
> Relavent contents can be found in Yade documentation:
> https://yade-dem.org/doc/user.html#scene-construction
> (user's manual--->scene construction-->creating interactions)
>
> I tried to run a simpler triaxial test using CohFrictMat as a first
> step. There are 500 particles, and parameters goes as follows: young =
> 5e6, normalcohesion = shearcohesion = 4500, confining pressure = 10kPa.
> Particles are packed randomly in a 0.15*0.15*0.3 box. Confining pressure
> is reached by internal compaction.
>
> After a stable status is reached, count the number of all interactions.
> Then the aabbenlargefunctor and interactiondetection factor is altered to
> 1.25(or 1.5). Next step, count the total number of interactions again and
> the two factors are set to 1.
> This method didn't work correctly:
>
> 1. Whatever the two factor is, 1.25 or 1.5, total interactions before and
> after alteration show little difference.
> 2. Besides, after the alteration, the utils.avgNumInteractions() is
> printed instantly and always turn out to be less than 7. According to [1],
> this number should be around 10 when intRaius is 1.25, and around 14 when
> intRadius is 1.5.
> 3. Stress and strain data is recorded, stress-strain curve show no obvious
> difference under different intRadius.
>
> The code is like this:
>
> *********************************************************************
> from yade import pack
>
> num_intr = 0
> O.materials.append(CohFrictMat(young=5e6, poisson=0.333,frictionAngle=0,normalCohesion
> = 4500, shearCohesion = 4500, density=2600, label='spheres'))
> O.materials.append(FrictMat(young=5e6, poisson=0.333,frictionAngle=0,
> density=0,label='walls'))
> O.bodies.append(aabbWalls([Vector3(0,0,0),Vector3(0.15,0.
> 3,0.15)],thickness=0,material='walls'))
> sp = pack.SpherePack()
> sp.makeCloud(Vector3(0,0,0),Vector3(0.15,0.3,0.15),-1,0.3333,500)
> O.bodies.append([sphere(center,rad,material='spheres') for center,rad in
> sp])
> triax=TriaxialStressController(maxMultiplier=1.+4e-3,
> finalMaxMultiplier=1.+1e-4,stressMask = 7,internalCompaction=True)
> triax.goal1=triax.goal2=triax.goal3=-1e4
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(label =
> 'bo1s'),Bo1_Box_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom6D(label =
> 'ig2ss'),Ig2_Box_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label
> = 'ip2cc')],
> [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_
> CohFrictPhys_CohesionMoment()]
> ),
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,
> timestepSafetyCoefficient=0.8),
> triax,
> NewtonIntegrator(damping=0.2),
> PyRunner(command = 'isotropic_compaction()', iterPeriod = 1000, label
> = 'checker'),
> ]
>
> def isotropic_compaction():
> global num_intr
> if unbalancedForce()<0.01 and abs(-1e4-triax.meanStress)/1e4<0.001:
> for i in O.interactions:
> num_intr = num_intr + 1
> print '*************************************************'
> print 'total interactions before reset intR is ', num_intr
> print 'current step is ',O.iter
> bo1s.aabbEnlargeFactor = 1.25
> ig2ss.interactionDetectionFactor = 1.25
> ip2cc.setCohesionNow = True
> num_intr = 0
> checker.iterPeriod = 1
> checker.command = 'set_interaction()'
> O.pause()
>
> def set_interaction():
> global num_intr
> bo1s.aabbEnlargeFactor = -1
> ig2ss.interactionDetectionFactor = -1
> for i in O.interactions:
> num_intr = num_intr + 1
> print '************************************************'
> print 'current step is ',O.iter
> print 'coordn = ', utils.avgNumInteractions()
> print 'total interactions after reset intR is ', num_intr
> num_intr = 0
> O.pause()
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Thank Jan Stránský for carefully reading and testing the scipt. I was hoping to provide a script that can easily run, so the key parts relavent to intR may not seem emphasized. By "before reset intR" in the script, I actually meant "before setting intR to 1.25", sorry for the language mistake :D

To prevent confusion may I confirm this sentence:"So you will get actual number of interactions in set_interactions() function, but the number will be before reset of intR". The number of interactions I get in set_interaction() function is the number of interactions when intRadius is 1.25, is that right?

In fact I've also noticed the 4500e10 matter. Just as Jérôme Duriez pointed, I ran a small script, in which 6 spheres(CohFrictMat) are ranged like a 2*3 matrix. Two of them are bigger than others and are in contact with each other even when intR is 1.0. Other four spheres are smaller, do not touch any other spheres at the very beginning. So if intR is 1.0, only two big spheres are in contact, num_intr is 1("num_intr" is short version of"number of interactions", and "intR" is short version of "intRadius")

So I set the intR big enough for the smaller spheres to contact and count num_intr, it still remains 1. No new contacts are created. I tried by adding "setCohesionNow=True" when intR is altered, SOMETHING happens:
1. If normalCohesion of the material is a smaller value(like 4500), still nothing.
2. if normalCohesion is bigger(like 450000), two new contacts are created between the bigger sphere and the nearest small sphere. And the small spheres will be "stick" to the bigger ones.
3. If normalCohesion is even bigger(like 500000), five new contacts are created between all neighbouring spheres. And they are all "stick" together after setting cohesion.

I've got two questions about this phenomena:
1. Why is contact detection rely on mechanical behavior(In this case different normalCohesion value affect interaction number)?
2. In my understanding, setCohesionNow() merely set a maximum tensile strength between two particles(not a strength, but a maximum strength value), but in this occasion, after setting cohesion, some "strength" are generated between particles, why is that?
As a first guess, this phenomena may have something to do with the CohFrictMat I used. Is there an "expert" in this material who can explain this?

Forgot to paste the script. The testing script mentioned in #4 is like this:

**********************************
O.materials.append(CohFrictMat(young = 1e6, poisson = 0.5,normalCohesion = 500000, shearCohesion = 500000, frictionAngle = 0,density = 2700,label='spheres'))

O.bodies.append([
    utils.sphere(center=(0.01,0,0.01),radius=0.008,material = 'spheres'),
    utils.sphere(center=(-0.01,0,0.01),radius=0.008,material = 'spheres'),
    utils.sphere(center=(0.01,0,0.03),radius=0.008,material = 'spheres'),
    utils.sphere(center=(-0.01,0,0.03),radius=0.008,material = 'spheres'),
    utils.sphere(center=(0.01,0,-0.01),radius=0.01,fixed=True,material = 'spheres'),
    utils.sphere(center=(-0.01,0,-0.01),radius=0.01,fixed=True,material = 'spheres'),
])

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor = 1.5,label = 'bo1s')]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor = 1.5,label = 'ig2ss')],
        [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label = 'ip2cc')],
        [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
    ),
    NewtonIntegrator(damping=0.1),
    PyRunner(command = 'stick()',iterPeriod = 10000,label = 'sticker')
]

ip2cc.setCohesionNow = True
O.step()
num_intr = 0
print "current step is ", O.iter
for i in O.interactions:
    num_intr = num_intr + 1
print 'interaction number is ', num_intr
num_intr = 0

def stick():
    global num_intr
    print "current step is ", O.iter
    for i in O.interactions:
        num_intr = num_intr+1
    print "interaction number is ", num_intr
    num_intr = 0
    sticker. iterPeriod = 100
    O.pause()

O.dt=0.005*utils.PWaveTimeStep()
*****************************************

And some questions about "launchpad" functions:
1. How can I answer to a specific person(or “floor” like #1,#2, etc), like Jan Stránský did in #3?
2. Would it be possible to invite a specific person to answer my question? I ask this because Luc Scholtès is the author of paper [1] mention in my original question and I saw him in this "Yade answer forum" too. It would be great if he could offer some advice.

Thanks again to dear Jérôme Duriez and Jan Stránský for their precious help.

Hi,

> 1. Why is contact detection rely on mechanical behavior(In this case different normalCohesion value affect interaction number)?

If cohesion is small the interactions are created then immediately broken for the tensile force exceeds the tensile strength.
Thus one possible trick is to set high cohesion to create the ints, then set the relevant interaction->adhesion at the next step within a loop.
We could add one option to initialize the force-free distance as the initial distance, it would make this particular situation a bit simpler (I think it is done this way in another (CPM?) contact model).

> 2. In my understanding, setCohesionNow() merely set a maximum tensile strength between two particles(not a strength, but a maximum strength value), but in this occasion, after setting cohesion, some "strength" are generated between particles, why is that?

To me "strength" means "maximum force". Not sure what "maximum strength" is. If you mean "some forces are generated", that is correct: the distant interactions are initially in traction (regardless of the strength except if it exceeds strength, in which case it is null).

> 1. How can I answer to a specific person(or “floor” like #1,#2, etc), like Jan Stránský did in #3?

"Reply to" the question from an email client, without erasing the replied-to (more a "lazy" than a "good" practice actually).

> 2. Would it be possible to invite a specific person

IIRC lanchpad lets you send 5 messages/day to other lauchpad accounts (click on the launchpad name -> contact this user).

Bruno

>Thus one possible trick is to set high cohesion to create the ints, then set the relevant interaction->adhesion at the next step within a loop.

I forgot to add: make the initial configuration force-free if that's what you need, with i.phys.unp = i.geom.penetrationDepth.

Also, you may note that the aforementioned paper [1] used in fact "JCFpm" model (see [*] and related classes) to describe particles interactions (rather than CohFrictMat).

As a particular feature, this JCFpm model uses the very initial particle distance as the force free distance (corresponding to unp in CohFrictMat model)

If you're sticking for a while with the CohFrictMat model, note also YADE graphical user interface offers a sub-step option (just check the corresponding box) that may allow you to check if the distant interaction exists at some point within a time step (before being possibly immediately erased by the contact law in the same time step, as explained by Bruno in #7)

[*] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Ip2_JCFpmMat_JCFpmMat_JCFpmPhys

>"YADE graphical user interface offers a sub-step option (just check the corresponding box) that may allow you to check if the distant interaction exists at some point within a time step (before being possibly immediately erased by the contact law in the same time step"

I don't think it is possible because everything happens in interactionLoop (one single sub-step).

Thank Bruno Chareyre and Jerome!
When I said "strength" in #4 I actually meant "force", sorry for the language mistake.

I tried to use this JCMpm model. But there is still something wrong. Firstly the isotropic state is reached by internal compaction. Then the cohesiveTresholdIteration is set two steps ahead of current step, so that new contact could be set cohesive. Meanwhile intR is set to 1.5. But when cohesive contact is checked at next step, it turn out that:
1. Only new contacts(created by sphere enlarge factor) is set cohesive. Other contacts are still frictional. Is there a way to set all interaction between particles cohesive?
2. When intR is set to 1.5, the coordination number remains to be about 9.8(almost the same with intr=1.25). Expected coord number should be 13. How might this be explained?

More about question 1. I tried to set the cohesiveTresholdIteration a big number(like 100000) to ensure that this number exceeds the compaction-finish timestep, then reset it as soon as compaction process finish. But in this way, some cohesive contacts break during compaction process. So it might be best if all contacts can be set within one single timestep.

Here is the code. Many thanks!

*************************************************
from yade import pack

num_intr = 0

O.materials.append(JCFpmMat(young=5e6, poisson=0.333,frictionAngle=0,cohesion = 4500, tensileStrength = 4500, density=2600, label='spheres'))
O.materials.append(FrictMat(young=5e6, poisson=0.333,frictionAngle=0,density=0,label='walls'))

O.bodies.append(aabbWalls([Vector3(0,0,0),Vector3(0.15,0.3,0.15)],thickness=0,material='walls'))
sp = pack.SpherePack()
sp.makeCloud(Vector3(0,0,0),Vector3(0.15,0.3,0.15),-1,0.3333,500)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

triax=TriaxialStressController(maxMultiplier=1.+4e-3,finalMaxMultiplier=1.+1e-4,stressMask = 7,internalCompaction=True)
triax.goal1=triax.goal2=triax.goal3=-1e4

O.engines=[
    ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(label = 'bo1s'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(label = 'ig2ss'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(label = 'ip2jj')],
  [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    NewtonIntegrator(damping=0.2),
    PyRunner(command = 'isotropic_compaction()', iterPeriod = 1000, label = 'checker'),
]

#############################
#######Key Part Below############
#############################
def isotropic_compaction():
    global num_intr
    if unbalancedForce()<0.01 and abs(-1e4-triax.meanStress)/1e4<0.001:
        for i in O.interactions:
            if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
                num_intr += 1
        triax.internalCompaction=False
        print '*************************************************'
        print 'current step is ',O.iter
        print 'total interactions before reset intR is ', num_intr
        ip2jj.cohesiveTresholdIteration = O.iter + 2
        bo1s.aabbEnlargeFactor = 1.5
        ig2ss.interactionDetectionFactor = 1.5
        num_intr = 0
        checker.iterPeriod = 1
        checker.command = 'set_interaction()'
        O.pause()

def set_interaction():
    O.pause()
    print '************************************************'
    print 'current step is ',O.iter
    bo1s.aabbEnlargeFactor = -1
    ig2ss.interactionDetectionFactor = -1
    print 'coordn = ',utils.avgNumInteractions()
    num_intr=0
    num_coh=0
    for i in O.interactions:
        if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
            num_intr += 1
            if i.phys.isCohesive:
                num_coh +=1
    print 'num_coh = ', num_coh
    print 'total interactions after set intR is ', num_intr

1. The behavior you observed is normal. As the doc of cohesiveTresholdIteration [1] says, Ip2_JCFpmMat_JCFpmMat_JCFpmPhys will set as cohesive only the *new* contacts.

In fact, it is a general rule in YADE that any Ip2_* executes his job (which is to define interaction properties) only once, where interactions are created.

If you'd like to "apply" cohesion to a previously created packing, you have to perform yourself Ip2_JCFpmMat_JCFpmMat_JCFpmPhys's corresponding job, using Python loops:

for interaction in O.interactions:
   interaction.phys.isCohesive = True
   # define cohesive thresholds interaction.phys.FnMax [2] and interaction.phys.FsMax [3] from material properties such as done in source code [4]. See also corresponding equations e.g. in [Duriez2016] (listed among the references of the doc)

2. Why are you expecting 13 as a coordination number ? The coordination number you will eventually get for a given interaction radius highly depends on the packing you're considering.
There is no unique relationship coordination number = f(interaction radius)

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Ip2_JCFpmMat_JCFpmMat_JCFpmPhys.cohesiveTresholdIteration
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.JCFpmPhys.FnMax
[3 ]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.JCFpmPhys.FsMax
[4] https://github.com/yade/trunk/blob/master/pkg/dem/JointedCohesiveFrictionalPM.cpp#L260 and line below

Problem solved. Thank Jérôme Duriez. Now I know how to apply cohesion.

Also thank Bruno Chareyre and Jan Stránský :D