------------------------------------------------------------ revno: 1755 committer: Václav Šmilauer <vac...@flux> branch nick: trunk timestamp: Wed 2009-09-02 17:40:53 +0200 message: 1. Add scene that has been sitting here for a few days (should be run for comparison with pfc3d) added: examples/cyl.py
-- lp:yade https://code.launchpad.net/~yade-dev/yade/trunk Your team Yade developers is subscribed to branch lp:yade. To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.
=== added file 'examples/cyl.py' --- examples/cyl.py 1970-01-01 00:00:00 +0000 +++ examples/cyl.py 2009-09-02 15:40:53 +0000 @@ -0,0 +1,56 @@ +# encoding: utf-8 +cylHt,cylRd=1,.2 +nSpheres=2e4 + +def unitSquare(): + """Return square composed of 2 facets in the xy plane, centered at zero with unit extents.""" + import gts + vv=[gts.Vertex(1,1,0),gts.Vertex(-1,1,0),gts.Vertex(-1,-1,0),gts.Vertex(1,-1,0)] + ee=[gts.Edge(vv[0],vv[1]),gts.Edge(vv[1],vv[2]),gts.Edge(vv[2],vv[3]),gts.Edge(vv[3],vv[0]),gts.Edge(vv[0],vv[2])] + surf=gts.Surface() + surf.add(gts.Face(ee[0],ee[1],ee[4])); surf.add(gts.Face(ee[2],ee[3],ee[4])) + return surf + +def unitCylinder(nDiv=24): + """Returns GTS surface approximating cylinder (without caps), with + of height 2 and radius 2, centered at origin, axis coincident with + the z-axis. + + @param nDiv: polyhedron approximating circle. + """ + import numpy; from yade import pack + thetas=numpy.linspace(0,2*pi,nDiv,endpoint=True) + ptsBase=[Vector3(cos(th),sin(th),-1) for th in thetas] + ptsTop=[p+Vector3(0,0,2) for p in ptsBase] + return pack.sweptPolylines2gtsSurface([ptsBase,ptsTop]) + +from yade import pack,timing +cyl=unitCylinder(); sq=unitSquare(); sq.translate(0,0,-1); cyl.copy(sq) +cyl.scale(cylRd,cylRd,.5*cylHt); cyl.rotate(1,0,0,-pi/4) # 45° anti-colckwise in the yz plane +# calling gtsSurface2Facets with just "cyl" (without constructing the faces tuple) ignores 2 faces that were copy'd before; bug in pygts? +cylIds=O.bodies.append(pack.gtsSurface2Facets(cyl.faces())) +sp=pack.SpherePack(); wd=cylRd*sqrt(2); rMean=(.2*wd*wd*cylHt/(nSpheres*(4/3.)*pi))**(1/3.) +print 'Generating cloudâ¦' +sp.makeCloud((-wd/2,-wd/2,-.5*cylHt),(wd/2,wd/2,.5*cylHt),rMean,0,int(nSpheres),False) +sp.rotate((1,0,0),-pi/4) +O.bodies.append([utils.sphere(s[0],s[1],density=3000) for s in sp]) + +O.engines=[ + BexResetter(), + BoundingVolumeMetaEngine([InteractingSphere2AABB(),InteractingFacet2AABB(),MetaInteractingGeometry2AABB()]), + InsertionSortCollider(nBins=5,sweepLength=.1*rMean), + InteractionDispatchers( + [ef2_Sphere_Sphere_Dem3DofGeom(),ef2_Facet_Sphere_Dem3DofGeom()], + [SimpleElasticRelationships()], + [Law2_Dem3Dof_Elastic_Elastic()], + ), + GravityEngine(gravity=(0,0,-1e3)), # gravity artificially high, to make it faster going ;-) + RotationEngine(rotateAroundZero=True,zeroPoint=(0,0,0),rotationAxis=(0,1,1),angularVelocity=30*(2*pi/60),subscribedBodies=cylIds,label='rotor'), + NewtonsDampedLaw(damping=.3), +] +O.dt=utils.PWaveTimeStep() +O.stopAtIter=int(2*pi/(rotor['angularVelocity']*O.dt)) +O.timingEnabled=True; timing.reset() + + +
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
