yade-users team mailing list archive

[nobody <question690079@xxxxxxxxxxxxxxxxxxxxx>]

New question #690079 on Yade:
https://answers.launchpad.net/yade/+question/690079

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

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.