← Back to team overview

yade-users team mailing list archive

Re: [Question #706962]: The problem of matchmaker

 

Question #706962 on Yade changed:
https://answers.launchpad.net/yade/+question/706962

William posted a new comment:
Hi, Jan. I have simplified the code.

Yade 2021.01a

> Here is the new code and the error.

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)

## 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
sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=AvgRadius,rRelFuzz=0,num=num_spheres,distributeMass=False,seed=1)
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)
    b.state.mass*=mat1.density/mat1.density
    sphere1.append(b.id)
  else:
    b.mat = mat2
    b.shape.color = (0,1,1)
    b.state.mass*=mat2.density/mat1.density
    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=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
]


#######################################
###   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      ###")
  print('----------Time:'+ str((time.time()-time0)/60) + 'min----------')

> error:

terminate called after throwing an instance of 'std::invalid_argument'
  what():  MatchMaker: no match for (0,1), and values required for algo computation 'avg' not specified.

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