problem with createInteraction() command

Asked by Son Pham Thai on 2018-12-05

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()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bruno Chareyre
Solved:
2018-12-05
Last query:
2018-12-05
Last reply:
2018-12-05
Jérôme Duriez (jduriez) said : #1

Hi,

Computers are always right: if YADE says "Interaction already exists", the interaction most certainly exists.

A possible scenario explaining why you can not detect this interaction could be that this interaction is not real yet (b.intrs() gives real interactions only, for instance).

 Could you check what

 O.interactions.has(39,18)

says ? (this would be True for a not real interaction as well)

Hi Son,
As Jerome suggests the apparent contradiction is because the interaction exists in a virtual state: detected by collider, not activated by contact model. It is more likely to happen with enlarged bounds ,but it can happen in the general case to.
How to create an interaction when it is already in a virtual state is actually a good question. I'm not sure how to solve it yet.

I'll check if it can be fixed in the code.

Bruno

Hi Jérôme,

Thank you for your response.

O.interactions[39,18] returns this: <Interaction instance at 0x7f20cc04a5b0>

I am reading about this real and not-real-yet interaction in Yade documentation now.

Cheers,
Thai-Son

Thanks Bruno Chareyre, that solved my question.

Son,
You script runs smoothly after a change [1] in the source code.
yadedaily should be updated soon.
Bruno

[1] https://github.com/yade/trunk/commit/20f4cf79ed7c4a92659e423ddbc775af7ee1394d

Hi Bruno,

Thank you very much for the new commit ;)

Cheers,
Son.