Create a pack in which the cylindrical particles and spherical particles are randomly distributed.

Hello,

a quick solution would be:
1) create ordinary sphere packing using makeCloud (possible called multiple times with different parameters for spheres and spheres to be replaced by cylinders)
2) replace respective spheres by cylinders
3) (if needed) compress the packing

For more detailed answer, please provide more detailed description of the problem:
- if/what you have already tried
- definition of "random distribution" (positions? orientation? volume fraction? ...)
- ...

cheers
Jan

Thank you very much for your prompt reply, Jan.

1. For your questions:

(1) I already created the cylindrical particles. My code is as follows:

##############################################

# encoding: utf-8
"An example showing how to create cylinders with random length."

from yade.gridpfacet import *
from yade import pack, plot,qt
import gts, os.path, locale

locale.setlocale(locale.LC_ALL, 'en_US.UTF-8') # Note: gts is locale-dependent. If, for example, german locale is used, gts.read()-function does not import floats normally

#### Parameters ####
Lmin=3. # minimum length
Lmax=6. # maximum length
n=10 # number of cylinders in one direction
r=0.5 # radius of the cylinder element
phi=30. # friction angle
E=1e6 # Young's modulus
radius=1e-02
wire=False
fixed=False
z=-1.2
color=[0,0,1]
Ns=200
R=0.1 #radius of the cylinder
L=5 #length of the cylinder
dx=0.5
dy=0
nx=20
ny=1
z=2
frictionAngle1=radians(15)
frictionAngle2=radians(15)
frictionAngle3=radians(15)

#### Engines ####
O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_GridConnection_Aabb(),
  Bo1_PFacet_Aabb(),
  Bo1_Sphere_Aabb(),
 ]),
 InteractionLoop([
  Ig2_GridNode_GridNode_GridNodeGeom6D(),
  Ig2_Sphere_Sphere_ScGeom(),
  Ig2_Sphere_GridConnection_ScGridCoGeom(),
  Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), # cylinder-cylinder interaction
  Ig2_Sphere_PFacet_ScGridCoGeom(), # needed for GridNode-pFacet interaction (why is this not included in Ig2_GridConnection_PFacet_ScGeom???)
  Ig2_GridConnection_PFacet_ScGeom(), # Cylinder-pFcet interaction
 ],
 [
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),
  Ip2_FrictMat_FrictMat_FrictPhys()
 ],
 [
  Law2_ScGeom_FrictPhys_CundallStrack(), # contact law for sphere-sphere
  Law2_ScGeom6D_CohFrictPhys_CohesionMoment(), # contact law for "internal" cylider forces
  Law2_ScGridCoGeom_FrictPhys_CundallStrack(), # contact law for Cylinder-pFacet interaction
  Law2_GridCoGridCoGeom_FrictPhys_CundallStrack() # contact law for cylinder-cylinder interaction
 ]
 ),
 GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.1,label='ts'),
 NewtonIntegrator(gravity=(0.,0,-10),damping=0.5,label='newton'),
 RotationEngine(angularVelocity=1,rotationAxis=(0,1,0),ids=[0,1,2])+TranslationEngine(translationAxis=[1,0,0],velocity=0.5,ids=[0,1,2]),
 #VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=5000),
 #PyRunner(command='addPlotData()',iterPeriod=100),

]
O.dt=0.1*PWaveTimeStep()

#### Creat materials ####
O.materials.append( CohFrictMat( young=E,poisson=0.3,density=1000,frictionAngle=radians(phi),normalCohesion=1e10,shearCohesion=1e10,momentRotationLaw=True,label='cMat' ) ) #material to create the gridConnections
O.materials.append( FrictMat( young=E,poisson=0.3,density=1000,frictionAngle=radians(phi),label='fMat' ) ) # material for general interactions

#cylindrical roller: the IDs of cylinder must be [0, 1, 2] which are used for translation and rotation engines.
nodesIds=[]
cylIds=[]
cylinder((-8,-8,1.5),(-8,8,1.5),radius=0.5,nodesIds=nodesIds,cylIds=cylIds,color=[1,0,0],fixed=False,intMaterial='cMat',extMaterial='fMat')

#### Spheres ####
O.materials.append(FrictMat(young=4.0e6,poisson=0.5,frictionAngle=frictionAngle1,density=1600,label='spheremat'))
sp=pack.SpherePack()
Ns=sp.makeCloud(Vector3(-5,-5,3.0),Vector3(+5,+5,+6.0),0.3,0.2,Ns,False,0.8)
O.bodies.append([sphere(center,rad,material='spheremat') for center,rad in sp])

nodesIds1=[]
cylIds1=[]
#### Create cylinders ####
for i in range(0,nx):
 x=-5+i*dx
 for j in range(0,ny):
  y=-2.5+j*dy
  cylinder((x,y+L,z),(x,y,z),radius=R,nodesIds=nodesIds1,cylIds=cylIds1,
      color=[1,0,0],fixed=False,intMaterial='cMat',extMaterial='fMat')

# collect history of data which will be plotted
def addPlotData():
 # each item is given a names, by which it can be the unsed in plot.plots
 # the **O.energy converts dictionary-like O.energy to plot.addData arguments
 plot.addData(i=O.iter,t=O.time,)
 plot.saveDataTxt('bbb.txt.bz2')

####create Pfacet from gts file####
#O.materials.append(CohFrictMat(young=1e8,poisson=0.3,density=2650,frictionAngle=radians(20),normalCohesion=1e100,shearCohesion=1e100,momentRotationLaw=True,label='gridNodeMat'))
#O.materials.append(FrictMat(young=1e8,poisson=0.3,density=2650,frictionAngle=radians(20),label='pFacetMat'))
#nodesIds1,cylIds1,pfIds1 = gtsPFacet('box.gts',shift=(3.,0.,1.2),scale=1,radius=radius,wire=wire,fixed=False,materialNodes='gridNodeMat',material='pFacetMat',color=color)

#### Creat ground with pFacets ####
color=[255./255.,102./255.,0./255.]
n0=O.bodies.append( gridNode([-10,-10,0],r,wire=False,fixed=True,material='cMat',color=color) )
n1=O.bodies.append( gridNode([10,-10,0],r,wire=False,fixed=True,material='cMat',color=color) )
n2=O.bodies.append( gridNode([10,10,0],r,wire=False,fixed=True,material='cMat',color=color) )
n3=O.bodies.append( gridNode([-10,10,0],r,wire=False,fixed=True,material='cMat',color=color) )
O.bodies.append( gridConnection(n0,n1,r,color=color,material='fMat') )
O.bodies.append( gridConnection(n1,n2,r,color=color,material='fMat') )
O.bodies.append( gridConnection(n2,n0,r,color=color,material='fMat') )
O.bodies.append( gridConnection(n2,n3,r,color=color,material='fMat') )
O.bodies.append( gridConnection(n3,n0,r,color=color,material='fMat') )
O.bodies.append( pfacet(n0,n1,n2,wire=False,material='fMat',color=color) )
O.bodies.append( pfacet(n0,n2,n3,wire=False,material='fMat',color=color) )

#### For viewing ####
qt.View()
Gl1_Sphere.stripes=True

#### Allows to reload the simulation ####
O.saveTmp()

#####################################################

(2). "random distribution" means the positions and orientations.

2. For your solutions:

Because my spherical particles and cylindrical particles have total different shapes, how to replace the spheres with cylinders?

Thank you a lot!

a MWE:
#########
def randomOrientation(): # most likely not very well random, but for this example it is ok
   from random import random as r
   return Quaternion((r(),r(),r()),2*pi*r())

# create loose spherical packing using makeCloud
sp = yade.pack.SpherePack()
mi,ma = (0,0,0),(20,20,20)
nCyls,nSphs = 40,30
sp.makeCloud(mi,ma,rMean=2,rRelFuzz=.5,num=nCyls) # makeCloud for cylinders
sp.makeCloud(mi,ma,rMean=1,rRelFuzz=.5,num=nSphs) # makeCloud for spheres
for i,(pos,radius) in enumerate(sp):
   if i < nCyls: # add cylinder
      ori = randomOrientation()
      extents = (radius,.1*radius,.1*radius)
      O.bodies.append(box(pos,extents,ori))
   else: # add sphere
      O.bodies.append(sphere(pos,radius))
##########
there is no cylinder command in my yade, so I used box instead, but the idea should be clear.
cheers
Jan

Thank you very much for your code, Jan. That is exactly what I want.