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>