yade-users team mailing list archive

[Son Pham Thai <question676609@xxxxxxxxxxxxxxxxxxxxx>]

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

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


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