local coordination number

Asked by ytang

I want to get the local coordination number/mechanical coordination number for a local region. The following is the code, I print all the necessary information.
but the results obtained by this function: yade.utils.avgNumInteractions #### is different from the coordination number that I get.

#############################################
from yade import ymport
from yade import pack,export,geom
import itertools
from numpy import *
import numpy as np
from yade import pack, plot, export, utils
import math
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
O.bodies.append(ymport.text("fenqu.txt"))
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()]),
    PyRunner(command='checkUnbalanced()',realPeriod=2),
    #PyRunner(command='addPlotData()',iterPeriod=100),
    PyRunner(command='subbox()',iterPeriod=100),
    #PyRunner(command='stress_export()',iterPeriod=100),
    NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
]
O.dt=.5*PWaveTimeStep()
print(len(O.bodies))
def checkUnbalanced():
    if unbalancedForce()<.00001:
        O.pause()
#O.run(100000000,True)

########################## CN ###########################
def subbox():
    global ball_list1
    ball_list1 =[]
    global zero_contact_p
    zero_contact_p = 0
    global one_contact_p
    one_contact_p = 0
    global intrs
    intrs = 0
    global intrs_number
    intrs_number = 0
    coodN = avgNumInteractions(cutoff=0.0, skipFree=False, considerClumps=False)
    coodN1 = avgNumInteractions(cutoff=0.0, skipFree=True, considerClumps=False)
    for b in O.bodies:
        if 0< b.state.pos[2] <= 1:
            if isinstance(b.shape,Sphere):
                b.shape.color =Vector3(00,00,255) # blue
                m = b.id
                ball_list1.append(m)
                for bb in ball_list1:
                    intrs_number =len(O.bodies[bb].intrs())
                    intrs = intrs+intrs_number
                    if O.bodies[bb].intrs() == 0:
                        zero_contact_p = zero_contact_p+1
                    if O.bodies[bb].intrs() == 1:
                        one_contact_p = one_contact_p+1
                    particle_number = len(ball_list1)
                    anverage_coord = 2*intrs/particle_number
                    mechanical_average_coord = (2*intrs-one_contact_p)/(particle_number-zero_contact_p-one_contact_p)
    print ('coodN is:',coodN)
    print('coodN1 is:',coodN1)
    print('zero contact numbers:',zero_contact_p)
    print('one contact numbers:',one_contact_p)
    print('this is partilce numbers:',particle_number)
    print('ball-list1:',ball_list1)
    print("this is the number in region 1:",len(ball_list1))
    print('this is intrs_number:',intrs_number)
    print('the total number of inters is:',intrs)
    print('this is the coordination number:',anverage_coord)
    print('this is the mechanical coordination number:',mechanical_average_coord)

the results:
('coodN is:', 2.125)
('coodN1 is:', 3.3333333333333335)
('zero contact numbers:', 0)
('one contact numbers:', 0)
('this is partilce numbers:', 6)
('ball-list1:', [10, 11, 12, 13, 14, 15])
('this is the number in region 1:', 6)
('this is intrs_number:', 4)
('the total number of inters is:', 79) ###### it seems the total interaction numbers is not correct. but I didn't find the error.
('this is the coordination number:', 26)
('this is the mechanical coordination number:', 26)

thanks!
Yong

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
ytang (ytang116) said :
#1

sorry for the 'fenqu.txt'
here is the former part code:
from yade import pack, plot, export
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=.2)
sp.toSimulation()
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()]),
    PyRunner(command='checkUnbalanced()',realPeriod=2),
    #PyRunner(command='addPlotData()',iterPeriod=100),
    PyRunner(command='subbox()',iterPeriod=10000),
    #PyRunner(command='stress_export()',iterPeriod=100),
    NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
]
O.dt=.5*PWaveTimeStep()
print(len(O.bodies))
def checkUnbalanced():
    if unbalancedForce()<.01:
        O.pause()

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

Hello,

> for b in O.bodies:
> ...
> ball_list1.append(m)
> for bb in ball_list1:

what is the intention of this? Shouldn't it be

for b in O.bodies:
   ...
      ball_list1.append(m)
for bb in ball_list1:

i.e. first collect relevant bodies in ball_list1 and THEN loop over them?

cheers
Jan

Revision history for this message
ytang (ytang116) said :
#3

Hi Jan,

what I want to is the same as you said, I want to get the coordination number for a local region, (as well as the local fabric tensor, macroscale stresses, so that I can bridge the microscale contact force with the macroscale stresses. for instance, P-Q plane in geomechanics. )

my first step is to collect the spheres in that region and then loop over that region to get the relevant microscale variables. ## it seems that in my code, I don't need to times 2 for the interactions, because when I loop for all the particles, the interactions already counted twice.

sorry for my bad python programming skills.

I have another question.
why the YADE didn't develop the class that can calculate the local variables, for instance, local coordination number, local fabric tensor, just like the measurement sphere function in PFC.

recently, I'm working on this, but my programming skill is poor.

best,
Yong

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

> sorry for my bad python programming skills.

no problem, everybody begins somehow :-)

> why the YADE didn't develop the class that can calculate the local variables, for instance, local coordination number, local fabric tensor, just like the measurement sphere function in PFC.

simply: nobody did need it / find it useful enough to make it a library function

There are some, e.g. utils.bodyStressTensors..

Filtering "local" bodies can be done easily using predicates:
###
predicate = pack.inAlignedBox(...) # or other instance or combination of several predicates
localBodies = [b for b in O.bodies if predicate(b.state.pos,b.shape.radius)]
###

The same for interactions (e.g. based on their contactPoint):
###
predicate = pack.inAlignedBox(...)
localIntrs = [i for i in O.interactions if predicate(i.geom.contactPoint)]
###

Computing coordination number, fabric tensor and other quantities (for each local body of from local interactions) is a few lines of code.

Averaging is just sum/amount..

So IMO the power of Python and already existing Yade functionality is the reason why Yade does not have these functions "builtin"..

> recently, I'm working on this, but my programming skill is poor.

it this is meant that you would like to add those functions to Yade, of course feel free. The fact that they are not present currently does not mean they would not be useful currently or in the future.

As an example, subbox "pythonized" (e.g. not to lost in the for loops).
###
def subbox():
   ball_list1 = [b for b in O.bodies if 0< b.state.pos[2]<= 1 and isinstance(b.shape,Sphere)]
   bintrs = [b.intrs() for b in ball_list1]
   intrs = sum(len(ii) for ii in bintrs) # (?!?!?!?)
   zero_contact_p = len([ii for ii in bintrs if len(intrs) == 0])
   one_contact_p = len([ii for ii in bintrs if len(intrs) == 1])
   particle_number = len(ball_list1)
   mechanical_average_coord = ...
###

Have a look at the interaction computation. For each body, you count some interactions twice (if both interaction bodies are contained in ball_list1), but some not (those "from local to outside").

cheers
Jan

Revision history for this message
ytang (ytang116) said :
#5

Hi Jan,

thank you very much!

your example helps me a lot. which means I don't need to use the sphere's position to set the local area. I can use the predicate function. that's nice.

best,
Yong

Revision history for this message
ytang (ytang116) said :
#6

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