problem with createInteraction() command

Hi all,

When I used createInteraction() command to manually create intrs between 2 spheres, I got an error:

RuntimeError: Interaction #39+#18 already exists.

I have then check these 2 spheres #39 and #18:

Sphere #39 has no intrs, len(O.bodies[39].intrs()) = 0, while sphere #18 has 7 intrs but #39 is not in the neighboring list of sphere #18.

Please see the MWE (Minimal Working Example™) generating this error down below. Please let me know why I have this RuntimeError error?

Best regards,
Thai-Son

PS: I have also tried aabbEnlargeFactor and interactionDetectionFactor but I could not obtain what I want, that's why I have to try doing this createInteraction manually. To be more specified, what I want is that all spheres in the packing has at least 1 intrs.

#-------------- MWE ------------------------------#
# encoding: utf-8
# ------------------------------------------------------------------------------ #

#######################
### INPUT PARAMETERs ###
#######################

seed0=1 # = 1 -> same packing everytime
distributeMass0=False # True -> Mass fraction, False -> PARTICLE COUNT fraction

num_spheres0=200
mx=Vector3(.4e-2,.6e-2,.4e-2) # size of aggregate
mn=Vector3(0,0,0)
porosity0=0.4
rRelFuzz0=0.7 # for the case of polydisperse

compFricDegree = 15.0 # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
confiningS=-1e5
surfaceTension=0.07274 # [N/m]
young0= 50e9 # 50 GPa
poisson0= 0.2
density0 = 2530 # 2.53 g/cm3 = 2530 kg/m3
tensileStrength = 3.3e9 # 3.3 GPa
shearStrength = 3.3e9 # 3.3 GPa

####################################
## GENERATING PARTICLE PACKING ###
####################################
from yade import pack
import pylab

O.materials.append(CohFrictMat(young=young0,poisson=poisson0,density=density0,frictionAngle=radians(compFricDegree), normalCohesion=tensileStrength, shearCohesion=shearStrength, momentRotationLaw=True,etaRoll=0.1,label='spheres'))
O.materials.append(FrictMat(young=young0,poisson=poisson0,frictionAngle=0,density=0,label='frictionlessWalls'))

walls=aabbWalls((mn,mx),thickness=0,material='frictionlessWalls')
wallIds=O.bodies.append(walls)

### Aggregate generation ###

num_spheres1D=int(round(num_spheres0**(1./3))) # for calculation of rMean0
sp=pack.SpherePack()
rMean0=(1-rRelFuzz0) * ((mx[0]/num_spheres1D)/2) # Irregular case, in order to obtain: num_spheres_real not less than num_spheres
sp.makeCloud(minCorner=mn,maxCorner=mx,rMean=rMean0,rRelFuzz=rRelFuzz0,porosity=porosity0,num=num_spheres0,distributeMass=distributeMass0,seed=seed0)

O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

### mechanical loading #### To confer a solid-like nature to the packing

triax=TriaxialStressController(
 internalCompaction=True,
 goal1=confiningS,
 goal2=confiningS,
 goal3=confiningS,
 max_vel=10,
 label="triax"
)

newton=NewtonIntegrator(damping=0.4)

unsat=UnsaturatedEngine(dead=1,label="unsat")

# ---------------------------------------------------------------------------
O.engines=[
 ForceResetter(),
 ### Collsion detector
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),

 ### Interaction handling
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")],

  [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
   useIncrementalForm=True, #useIncrementalForm is turned on as we want plasticity on the contact moments
   always_use_moment_law=False, #if we want "rolling" friction even if the contact is not cohesive (or cohesion is broken), we will have to turn this true somewhere
   label='cohesiveLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.5),
 triax,
 newton
]

# ---------------------------------------------------------------------------
numDEMiter1stDEMloop=3000
numDEMiter2ndDEMloop=2000
numDEMiter3rdDEMloop=500
numDEMiterPerDryIter=250

while 1:
 O.run(numDEMiter1stDEMloop,True)
 unb=unbalancedForce()
 if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
  break
## REACH NEW EQU. STATE ###
finalFricDegree = 30 # contact friction during the deviatoric loading

# We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False
# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

while 1:
 O.run(numDEMiter2ndDEMloop,True)
 unb=unbalancedForce()
 if unb<0.001 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
  break

triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width

O.run(numDEMiter3rdDEMloop,True)
ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2]
si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2]

### counting the real number of spheres generated by YADE
num_spheres_real=0
for b in O.bodies:
 if not isinstance(b.shape,Box):
  num_spheres_real =num_spheres_real+1

num_spheres=num_spheres_real

### --------------------------------------------------------------------###
### ---- GLUEING THE PARTICLE manually ---- ###
### --------------------------------------------------------------------###

### calculating number of no-intrs bodies to see the effect of manuallly creating bond between particles

NoIntrsBodies=[]
for i in range(6,len(O.bodies)): # exclude walls, id<6
 if len(O.bodies[i].intrs()) == 0 and not i in NoIntrsBodies :
  NoIntrsBodies+=[i]

print 'NoIntrsBodies before creating additional bonds: ', len(NoIntrsBodies)

#-----------------------------------------------------------#
# CREATE BONDS BETWEEN PARTICLES MANUALLY
#-----------------------------------------------------------#
# List of no-bond particles

EnlargeFactorManual=150./100.

ParticleWithManualBondCreatedList=[]
for i in NoIntrsBodies:
 if not i in ParticleWithManualBondCreatedList:
  for j in range(6,len(O.bodies)): # exclude walls, id<6, only manually create bonds between spheres
   if not i==j and not i in ParticleWithManualBondCreatedList: # skip the sphere i itself and skip i if i has been already added in this "for" loop
    if (float((O.bodies[i].state.pos-O.bodies[j].state.pos).norm())-O.bodies[j].shape.radius) / O.bodies[i].shape.radius < EnlargeFactorManual:
    # --- MANUALLY CREATE INTERATION (BOND) BETWEEN 2 SPHERES --- #
     createInteraction(i,j)

     # -- this ParticleWithManualBondCreatedList is used for avoiding checking a sphere which already turned from no-bond to bonded sphere
     if not i in ParticleWithManualBondCreatedList:
      ParticleWithManualBondCreatedList+=[i]
     if not j in ParticleWithManualBondCreatedList:
      ParticleWithManualBondCreatedList+=[j]

NoIntrsBodiesAfter=[]
for i in range(6,len(O.bodies)): # exclude walls, id<6
 if len(O.bodies[i].intrs()) == 0 and not i in NoIntrsBodiesAfter :
  NoIntrsBodiesAfter+=[i]
print 'EnlargeFactorManual: ', EnlargeFactorManual
print 'NoIntrsBodies after creating additional bonds: ', len(NoIntrsBodiesAfter)

### --- DEBUGGING section --- # Please use this pack of codes after getting the error "RuntimeError: Interaction #39+#18 already exists."
#%cpaste
print ''
print "i in NoIntrsBodies >>>", i in NoIntrsBodies
print "j in NoIntrsBodies >>>", j in NoIntrsBodies
print "i in ParticleWithManualBondCreatedList >>>",i in ParticleWithManualBondCreatedList
print "j in ParticleWithManualBondCreatedList >>>",j in ParticleWithManualBondCreatedList
print "EnlargeFactorManual: ",EnlargeFactorManual
print "len(O.bodies[i].intrs()):", len(O.bodies[i].intrs())
print "len(O.bodies[j].intrs()):", len(O.bodies[j].intrs())
if len(O.bodies[i].intrs()) >0:
 for k in range(len(O.bodies[i].intrs())):
  print "O.bodies[i].intrs()[",k,"].dict(): id1=", O.bodies[i].intrs()[k].id1, " id2=", O.bodies[i].intrs()[k].id2
  print " ", O.bodies[i].intrs()[k].dict()
if len(O.bodies[j].intrs()) >0:
 for k in range(len(O.bodies[j].intrs())):
  print "O.bodies[j].intrs()[",k,"].dict(): id1=", O.bodies[j].intrs()[k].id1, " id2=", O.bodies[j].intrs()[k].id2
  print " ", O.bodies[j].intrs()[k].dict()