# yade-users team mailing list archive

## Re: [Question #670078]: Dense packing with given distribution of radii

```Question #670078 on Yade changed:

Jan Stránský proposed the following answer:
Hi,

> But I am unsure how I can make my spherePack dense after filling the predicate?
> And since my body of spheres is also spherical, it seems to me that a compression from 6 different planes will only distort my sphere.

one the approach of randomDensePack is to create a densePacking and than
"crop" it according to predicate

> it is suggested to change the source code for spherePack, how can this
be done?

e.g.:

####################################
def myMakeCloud(mn,mx,mu,sigma,rMin,num,maxTry=1000,seed=0): # [1]
import random
import numpy as np
random.seed(seed)
mn = Vector3(mn)
mx = Vector3(mx)
size = mx-mn
rnd = random.random
ret = []
for t in xrange(maxTry):
c = Vector3([mn[axis] + r + (size[axis]-2*r)*rnd() for axis in (0,1,2)])
packSize = len(ret)
overlap = False
for pc,pr in ret:
if pow(pr+r,2) >= (pc-c).squaredNorm():
overlap = True
break
if not overlap:
ret.append((c,r))
break
if t == maxTry:
raise RuntimeError, "WRNING: Exceeded {} tries to insert non-overlapping sphere to packing. Only {} spheres were added, although you requested {}.".format(maxTry,i,num)
return pack.SpherePack(ret)

def myRandomDensePack(predicate,mu,sigma,rMin,material=-1,spheresInCell=300,color=None,returnSpherePack=None,seed=0): # [2]
dim=predicate.dim()
if max(dim)==float('inf'): raise RuntimeError("Infinite predicate and no dimension of packing requested.")
center=predicate.center()
fullDim=dim
cloudPorosity=0.25
beta,gamma=fullDim[1]/fullDim[0],fullDim[2]/fullDim[0]
N100=spheresInCell/cloudPorosity
x1=mu*(1/(beta*gamma)*N100*(4/3.)*pi)**(1/3.)
y1,z1=beta*x1,gamma*x1; vol0=x1*y1*z1
O.switchScene(); O.resetThisScene() ### !!
O.periodic=True
O.cell.setBox((x1,y1,z1))
O.materials.append(FrictMat(young=3e10,density=2400))
sp = myMakeCloud(Vector3().Zero,O.cell.refSize,mu,sigma,rMin,spheresInCell,seed=seed)
O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=.05*mu),InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),PeriIsoCompressor(charLen=2*mu,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=1,keepProportions=True),NewtonIntegrator(damping=.6)]
O.materials.append(FrictMat(young=30e9,frictionAngle=.5,poisson=.3,density=1e3))
for s in sp: O.bodies.append(utils.sphere(s[0],s[1]))
O.dt=utils.PWaveTimeStep()
O.run(); O.wait()
sp=SpherePack(); sp.fromSimulation()
O.switchScene() ### !!
sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2]))
return filterSpherePack(predicate,sp,material=material,color=color,returnSpherePack=returnSpherePack)

mu, sigma = 20,2
rMin = .25*mu
sphs = myRandomDensePack(pred,mu,sigma,rMin,spheresInCell=300)
O.bodies.append(sphs)
####################################

cheers
Jan