View source for ScriptGenerateCylindricalPacking

<source lang="python">
from __future__ import division
from math import *
rCyl=1
hCyl=2
nPoly=20

hBox=hCyl/5 

phiStep=2*pi/nPoly

o=Omega()

from yade import utils
for n in range(nPoly):
	phi1,phi2=n*phiStep,(n+1)*phiStep
	def pt(angle,radius,z):
		return radius*sin(angle),radius*cos(angle),z
	a,b,c,d=pt(phi1,rCyl,0),pt(phi2,rCyl,0),pt(phi1,rCyl,hCyl),pt(phi2,rCyl,hCyl)
	o.bodies.append([
		utils.facet([a,b,c]),
		utils.facet([b,c,d])])
o.bodies.append(utils.box([0,0,-.5*hBox],[rCyl,rCyl,.5*hBox]))

for b in o.bodies: b['isDynamic']=False

rSphere=rCyl/10.
from numpy import arange
import random
for z in arange(rSphere,20*rSphere,2.1*rSphere):
	for r in arange(2*rSphere,rCyl-1.1*rSphere,2.1*rSphere):
		for theta in arange(0,2*pi-rSphere/r,2.2*rSphere/r):
			theta2=theta*(1+.01*(.5+random.random()))
			o.bodies.append(utils.sphere([r*sin(theta2),r*cos(theta2),z],rSphere*(1-.3*(random.random())),wire=True,color=[random.random(),random.random(),random.random()]))

o.initializers=[
	StandAloneEngine('PhysicalActionContainerInitializer'),
	MetaEngine('BoundingVolumeMetaEngine',[EngineUnit('InteractingSphere2AABB'),EngineUnit('InteractingFacet2AABB'),EngineUnit('MetaInteractingGeometry2AABB')])
	]
o.engines=[
	StandAloneEngine('PhysicalActionContainerReseter'),
	MetaEngine('BoundingVolumeMetaEngine',[
		EngineUnit('InteractingSphere2AABB'),
		EngineUnit('InteractingFacet2AABB'),
		EngineUnit('InteractingBox2AABB'),
		EngineUnit('MetaInteractingGeometry2AABB')
	]),
	StandAloneEngine('PersistentSAPCollider'),
	MetaEngine('InteractionGeometryMetaEngine',[
		EngineUnit('InteractingSphere2InteractingSphere4SpheresContactGeometry'),
		EngineUnit('InteractingFacet2InteractingSphere4SpheresContactGeometry'),
		EngineUnit('InteractingBox2InteractingSphere4SpheresContactGeometry')
	]),
	MetaEngine('InteractionPhysicsMetaEngine',[EngineUnit('SimpleElasticRelationships')]),
	StandAloneEngine('ElasticContactLaw'),
	DeusExMachina('GravityEngine',{'gravity':[0,0,-1e5]}),
	DeusExMachina('NewtonsDampedLaw',{'damping':0.4}),
	#StandAloneEngine('PeriodicPythonRunner',{'command':'uf=utils.unbalancedForce(); print uf;\nif uf<abs(grav["gravity"][2])*1e-9: o.pause()','realLim':2})

]
grav=[e for e in o.engines if e.name=='GravityEngine'][0]
o.dt=1e-6
o.stopAtIter=10000 # get some better stability criterion here
o.run()
#
o.wait()
utils.spheresToFile('/tmp/a.spheres')

</source>

[[Image:cyl-specimen-generation.png]]
[[Category:Python]]