gts-random-pack.py 3.19 KB
Newer Older
Christian Jakob's avatar
Christian Jakob committed
1 2 3 4 5

""" CAUTION:
Running this script can take very long!
"""

6 7 8 9
from numpy import arange
from yade import pack
import pylab
# define the section shape as polygon in 2d; repeat first point at the end to close the polygon
10
poly=((1e-2,5e-2),(5e-2,2e-2),(7e-2,-2e-2),(1e-2,-5e-2),(1e-2,5e-2))
11
# show us the meridian shape
12
#pylab.plot(*zip(*poly)); pylab.xlim(xmin=0); pylab.grid(); pylab.title('Meridian of the revolution surface\n(close to continue)'); pylab.gca().set_aspect(aspect='equal',adjustable='box'); pylab.show()
13
# angles at which we want this polygon to appear
14
thetas=arange(0,pi/2,pi/24)
15 16 17 18 19 20 21 22 23 24
# create 3d points from the 2d ones, turning the 2d meridian around the +y axis
# for each angle, put the poly a little bit higher (+2e-3*theta);
# this is just to demonstrate that you can do whatever here as long as the resulting
# meridian has the same number of points
#
# There is origin (translation) and orientation arguments, allowing to transform all the 3d points once computed.
#
# without these transformation, it would look a little simpler:
# 	pts=pack.revolutionSurfaceMeridians([[(pt[0],pt[1]+2e-3*theta) for pt in poly] for theta in thetas],thetas
#
25
pts=pack.revolutionSurfaceMeridians([[(pt[0],pt[1]+1e-2*theta) for pt in poly] for theta in thetas],thetas,origin=Vector3(0,-.05,.1),orientation=Quaternion((1,1,0),pi/4))
26 27 28 29 30 31 32
# connect meridians to make surfaces
# caps will close it at the beginning and the end
# threshold will merge points closer than 1e-4; this is important: we want it to be closed for filling
surf=pack.sweptPolylines2gtsSurface(pts,capStart=True,capEnd=True,threshold=1e-4)
# add the surface as facets to the simulation, to make it visible
O.bodies.append(pack.gtsSurface2Facets(surf,color=(1,0,1)))
# now fill the inGtsSurface predicate constructed form the same surface with sphere packing generated by TriaxialTest
33
# with given radius and standard deviation (see documentation of pack.randomDensePack)
34 35 36 37 38
#
# The memoizeDb will save resulting packing into given file and next time, if you run with the same
# parameters (or parameters that can be scaled to the same one),
# it will load the packing instead of running the triaxial compaction again.
# Try running for the second time to see the speed difference!
39
memoizeDb='/tmp/gts-triax-packings.sqlite'
40 41 42
sp=SpherePack()
sp=pack.randomDensePack(pack.inGtsSurface(surf),radius=5e-3,rRelFuzz=1e-4,memoizeDb=memoizeDb,returnSpherePack=True)
sp.toSimulation()
43 44 45
# We could also fill the horse with triaxial packing, but have nice approximation, the triaxial would run terribly long,
# since horse discard most volume of its bounding box
# Here, we would use a very crude one, however
46
if 1:
47
	import gts
48
	horse=gts.read(open('horse.coarse.gts')) #; horse.scale(.25,.25,.25)
49
	O.bodies.append(pack.gtsSurface2Facets(horse))
50 51
	sp=pack.randomDensePack(pack.inGtsSurface(horse),radius=5e-3,memoizeDb=memoizeDb,returnSpherePack=True)
	sp.toSimulation()
52 53
	horse.translate(.07,0,0)
	O.bodies.append(pack.gtsSurface2Facets(horse))
54
	# specifying spheresInCell makes the packing periodic, with the given number of spheres, proportions being equal to that of the predicate
55 56
	sp=pack.randomDensePack(pack.inGtsSurface(horse),radius=1e-3,spheresInCell=2000,memoizeDb=memoizeDb,returnSpherePack=True)
	sp.toSimulation()