The problem of matchmaker

When I add the 'MatchMaker' in the second engine, the error 'terminate called recursively' appears. I didn't find what's wrong with the code. If i delete the relavant code of matchmaker, it will run without the error.

Ubuntu 21.10

> Here is the error:
terminate called recursively
terminate called recursively
terminate called recursively
terminate called after throwing an instance of 'std::invalid_argument'

>Here is my code. this code is divided into two parts. if you are going to run it, you have to wait a moment (about 5 min) for finishing the first part and the start the 2nd part. the error is in the second part.

from yade import pack, timing
#time calculation
startT = O.time
endT = O.time
timeSpent = endT - startT
time0 = time.time()

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

nRead=readParamsFromTable(
 num_spheres=2000,
 key='_triax_base_',
 unknownOk=True)

from yade.params import table
num_spheres=table.num_spheres
key=table.key

damp=0.7
stabilityThreshold=0.001

## Sphere material
Young=300e6
Density=2650
Poisson=0.3
compFricDegree = 0
finalFricDegree = 27
normalCohesion=100e4
shearCohesion=100e4
etaRoll=0.1
AvgRadius=0.0005

## wall material
WYoung=6000e6
WPoisson=0.5
WDensity=0
WFrictionAngle=0

# corners of the initial packing
size=0.05
mn,mx=Vector3(0,0,0),Vector3(size,size,size)

## create materials for spheres and plates
O.materials.append(CohFrictMat(young = WYoung , poisson = WPoisson , frictionAngle = radians(WFrictionAngle) , density = WDensity,label='walls'))

mat1=CohFrictMat(young = Young , poisson = Poisson , frictionAngle = radians(compFricDegree) ,normalCohesion=100,shearCohesion=100,isCohesive=True,alphaKr=0,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True, density = Density,label='spheres1')
mat_1=O.materials.append(mat1)

mat2=CohFrictMat(young = Young , poisson = Poisson , frictionAngle = radians(compFricDegree) ,normalCohesion=normalCohesion,shearCohesion=shearCohesion,isCohesive=True,alphaKr=0,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True, density = Density,label='spheres2')
mat_2=O.materials.append(mat2)

print(mat_1,mat_2)

# mat_1=O.materials.append(CohFrictMat(young = Young , poisson = Poisson , frictionAngle = radians(compFricDegree) ,normalCohesion=100,shearCohesion=100,isCohesive=True,alphaKr=0,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True, density = Density,label='spheres1'))

# mat_2=O.materials.append(CohFrictMat(young = Young , poisson = Poisson , frictionAngle = radians(compFricDegree) ,normalCohesion=normalCohesion,shearCohesion=shearCohesion,isCohesive=True,alphaKr=0,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True, density = Density,label='spheres2'))

# O.materials.append(CohFrictMat(young = WYoung , poisson = WPoisson , frictionAngle = radians(WFrictionAngle) , density = WDensity,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
#psdSizes,psdCumm=[0.00099,0.00107,0.00116,0.00143,0.00163,0.00172,0.00209,0.00251,0.00424],[0.0,0.04,0.10,0.3,0.5,0.6,0.8,0.95,1]
sp=pack.SpherePack()
#sp.makeCloud(mn,mx,num=num_spheres,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,seed=1)
sp.makeCloud(mn,mx,rMean=AvgRadius,rRelFuzz=0,num=num_spheres,distributeMass=False,seed=1)
#sp.toSimulation(material='spheres1')
sp.toSimulation()

sphere1=[]
sphere2=[]

for b in O.bodies:
  if not isinstance(b.shape,Sphere): # change material only on spheres
    continue
  if random.random() < 0.9:
    b.mat = mat1
    b.shape.color = (1,0,0)
    sphere1.append(b.id)
  else:
    b.mat = mat2
    b.shape.color = (0,1,1)
    sphere2.append(b.id)

############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
 maxMultiplier=1.001,
 finalMaxMultiplier=1.00001,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
              [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
                [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=False,setCohesionOnNewContacts=False,label="cohesiveIp")],
                [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True,label='cohesiveLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]

#######################################
### APPLYING CONFINING PRESSURE ###
#######################################

triax.goal1=triax.goal2=triax.goal3=-90000

if 1:
  while 1:
    O.run(1000, True)
    unb=unbalancedForce()
    print ('unbalanced force:',unb,' mean stress: ',triax.meanStress," porosity:",triax.porosity,)
    if unb<stabilityThreshold and abs(-90000-triax.meanStress)/90000<0.001:
      break

  O.save('1.yade.bz2')
  print("### Sample saved ###")

  ###################################################
  ### REACHING A SPECIFIED POROSITY PRECISELY ###
  ###################################################

  triax.goal1=triax.goal2=triax.goal3=-100000
  triax.internalCompaction=False
  setContactFriction(radians(finalFricDegree))
  O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionOnNewContacts=True)]#,normalCohesion=MatchMaker(matches=((mat_1,mat_1,100),(mat_1,mat_2,normalCohesion/2),(mat_2,mat_2,normalCohesion))),shearCohesion=MatchMaker(matches=((mat_1,mat_1,100), (mat_1,mat_2,shearCohesion/2),(mat_2,mat_2,shearCohesion))))],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True,label='cohesiveLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
  ]

  while 1:
    O.run(1000, True)
    unb=unbalancedForce()
    print ('unbalanced force:',unb,' mean stress: ',triax.meanStress," porosity:",triax.porosity,)
    if unb<0.000005 and abs(-100000-triax.meanStress)/100000<0.001:
      break

  O.save('2.yade.gz')
  O.save('2.yade.bz2')
  print("### Isotropic state saved ###")
  print('----------Time:'+ str((time.time()-time0)/60) + 'min----------')