/usr/lib/python2.7/dist-packages/yade/pack.py is in python-yade 2017.01a-8.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 | # encoding: utf-8
# 2009 © Václav Šmilauer <eudoxos@arcig.cz>
"""
Creating packings and filling volumes defined by boundary representation or constructive solid geometry.
For examples, see
* :ysrc:`scripts/test/gts-operators.py`
* :ysrc:`scripts/test/gts-random-pack-obb.py`
* :ysrc:`scripts/test/gts-random-pack.py`
* :ysrc:`scripts/test/pack-cloud.py`
* :ysrc:`scripts/test/pack-predicates.py`
* :ysrc:`examples/packs/packs.py`
* :ysrc:`examples/gts-horse/gts-horse.py`
* :ysrc:`examples/WireMatPM/wirepackings.py`
"""
import itertools,warnings
from numpy import arange
from math import sqrt
from yade import utils
from yade.wrapper import *
from minieigen import *
## compatibility hack for python 2.5 (21/8/2009)
## can be safely removed at some point
if 'product' not in dir(itertools):
def product(*args, **kwds):
"http://docs.python.org/library/itertools.html#itertools.product"
pools = map(tuple, args) * kwds.get('repeat', 1); result = [[]]
for pool in pools: result = [x+[y] for x in result for y in pool]
for prod in result: yield tuple(prod)
itertools.product=product
# for now skip the import, but try in inGtsSurface constructor again, to give error if we really use it
try:
import gts
except ImportError: pass
# make c++ predicates available in this module
noPredicate = False
try:
from _packPredicates import * ## imported in randomDensePack as well
except ImportError: pass; noPredicate = True
# import SpherePack
from _packSpheres import *
from _packObb import *
##
# extend _packSphere.SpherePack c++ class by this method
##
if not (noPredicate):
def SpherePack_toSimulation(self,rot=Matrix3.Identity,**kw):
"""Append spheres directly to the simulation. In addition calling :yref:`O.bodies.append<BodyContainer.append>`,
this method also appropriately sets periodic cell information of the simulation.
>>> from yade import pack; from math import *
>>> sp=pack.SpherePack()
Create random periodic packing with 20 spheres:
>>> sp.makeCloud((0,0,0),(5,5,5),rMean=.5,rRelFuzz=.5,periodic=True,num=20)
20
Virgin simulation is aperiodic:
>>> O.reset()
>>> O.periodic
False
Add generated packing to the simulation, rotated by 45° along +z
>>> sp.toSimulation(rot=Quaternion((0,0,1),pi/4),color=(0,0,1))
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]
Periodic properties are transferred to the simulation correctly, including rotation (this could be avoided by explicitly passing "hSize=O.cell.hSize" as an argument):
>>> O.periodic
True
>>> O.cell.refSize
Vector3(5,5,5)
"""
#The following 2 lines do not work, because of accuaracy
#>>> O.cell.hSize
#Matrix3(3.53553,-3.53553,0, 3.53553,3.53553,0, 0,0,5)
"""
The current state (even if rotated) is taken as mechanically undeformed, i.e. with identity transformation:
>>> O.cell.trsf
Matrix3(1,0,0, 0,1,0, 0,0,1)
:param Quaternion/Matrix3 rot: rotation of the packing, which will be applied on spheres and will be used to set :yref:`Cell.trsf` as well.
:param \*\*kw: passed to :yref:`yade.utils.sphere`
:return: list of body ids added (like :yref:`O.bodies.append<BodyContainer.append>`)
"""
if isinstance(rot,Quaternion): rot=rot.toRotationMatrix()
assert(isinstance(rot,Matrix3))
if self.isPeriodic: O.periodic=True
if self.cellSize!=Vector3.Zero and self.isPeriodic:
O.cell.hSize=rot*Matrix3(self.cellSize[0],0,0, 0,self.cellSize[1],0, 0,0,self.cellSize[2])
O.cell.trsf=Matrix3.Identity
if not self.hasClumps():
return O.bodies.append([utils.sphere(rot*c,r,**kw) for c,r in self])
else:
standalone,clumps=self.getClumps()
ids=O.bodies.append([utils.sphere(rot*c,r,**kw) for c,r in self]) # append all spheres first
clumpIds=[]
userColor='color' in kw
for clump in clumps:
clumpIds.append(O.bodies.clump([ids[i] for i in clump])) # clump spheres with given ids together, creating the clump object as well
# make all spheres within one clump a single color, unless color was specified by the user
if not userColor:
for i in clump[1:]: O.bodies[ids[i]].shape.color=O.bodies[ids[clump[0]]].shape.color
return ids+clumpIds
SpherePack.toSimulation=SpherePack_toSimulation
class inGtsSurface_py(Predicate):
"""This class was re-implemented in c++, but should stay here to serve as reference for implementing
Predicates in pure python code. C++ allows us to play dirty tricks in GTS which are not accessible
through pygts itself; the performance penalty of pygts comes from fact that if constructs and destructs
bb tree for the surface at every invocation of gts.Point().is_inside(). That is cached in the c++ code,
provided that the surface is not manipulated with during lifetime of the object (user's responsibility).
---
Predicate for GTS surfaces. Constructed using an already existing surfaces, which must be closed.
import gts
surf=gts.read(open('horse.gts'))
inGtsSurface(surf)
.. note::
Padding is optionally supported by testing 6 points along the axes in the pad distance. This
must be enabled in the ctor by saying doSlowPad=True. If it is not enabled and pad is not zero,
warning is issued.
"""
def __init__(self,surf,noPad=False):
# call base class ctor; necessary for virtual methods to work as expected.
# see comments in _packPredicates.cpp for struct PredicateWrap.
super(inGtsSurface,self).__init__()
if not surf.is_closed(): raise RuntimeError("Surface for inGtsSurface predicate must be closed.")
self.surf=surf
self.noPad=noPad
inf=float('inf')
mn,mx=[inf,inf,inf],[-inf,-inf,-inf]
for v in surf.vertices():
c=v.coords()
mn,mx=[min(mn[i],c[i]) for i in 0,1,2],[max(mx[i],c[i]) for i in 0,1,2]
self.mn,self.mx=tuple(mn),tuple(mx)
import gts
def aabb(self): return self.mn,self.mx
def __call__(self,_pt,pad=0.):
p=gts.Point(*_pt)
if self.noPad:
if pad!=0: warnings.warn("Padding disabled in ctor, using 0 instead.")
return p.is_inside(self.surf)
pp=[gts.Point(_pt[0]-pad,_pt[1],_pt[2]),gts.Point(_pt[0]+pad,_pt[1],_pt[2]),gts.Point(_pt[0],_pt[1]-pad,_pt[2]),gts.Point(_pt[0],_pt[1]+pad,_pt[2]),gts.Point(_pt[0],_pt[1],_pt[2]-pad),gts.Point(_pt[0],_pt[1],_pt[2]+pad)]
return p.is_inside(self.surf) and pp[0].is_inside(self.surf) and pp[1].is_inside(self.surf) and pp[2].is_inside(self.surf) and pp[3].is_inside(self.surf) and pp[4].is_inside(self.surf) and pp[5].is_inside(self.surf)
class inSpace(Predicate):
"""Predicate returning True for any points, with infinite bounding box."""
def __init__(self, _center=Vector3().Zero): self._center=_center
def aabb(self):
inf=float('inf'); return Vector3(-inf,-inf,-inf),Vector3(inf,inf,inf)
def center(self): return self._center
def dim(self):
inf=float('inf'); return Vector3(inf,inf,inf)
def __call__(self,pt,pad): return True
class inHalfSpace(Predicate):
"""Predicate returning True any points, with infinite bounding box."""
def __init__(self, _center=Vector3().Zero, _dir=Vector3(1,0,0)):
self._center = Vector3(_center)
self._dir = Vector3(_dir)
assert self._dir.norm() > 0., "Direction has to be nonzero vector"
self._dir.normalize()
self._d = -self._dir.dot(self._center)
def aabb(self):
d,c = self._dir, self._center
inf=float('inf')
min = Vector3(-inf,-inf,-inf)
max = Vector3(+inf,+inf,+inf)
for i in xrange(3):
j = (i+1)%3
k = (i+2)%3
if d[i]==0 and d[j]==0:
if d[k] > 0: min[k] = c[k]
else: max[k] = c[k]
return min,max
def center(self):
return self._center
def dim(self):
inf=float('inf')
return Vector3(inf,inf,inf)
def __call__(self,pt,pad):
v = self._dir.dot(pt) + self._d
return v > pad
class inConvexPolyhedron(Predicate):
def __init__(self,planes):
self._inHalfSpaces = [inHalfSpace(c,d) for c,d in planes]
self._min, self._max = self._computeAabb()
def _computeAabb(self):
try:
import scipy.optimize
except ImportError:
raise ImportError, "scipy (package python-scipy) needed for pack.inConvexPolyhedron"
min,max = Vector3.Zero, Vector3.Zero
A,b = [],[]
for h in self._inHalfSpaces:
A.append(tuple(-h._dir))
b.append(h._d)
inf = float('inf')
for i in xrange(3):
c = Vector3.Zero
#
c[i] = 1
opt = scipy.optimize.linprog(c,A_ub=A,b_ub=b,bounds=(-inf,inf))
errmsg = "Something wrong in pack.inConvexPolyhedron defining planes.\nThe scipy.optimize.linprog output:\n{}\n"
if not opt.success:
raise ValueError, errmsg.format(opt)
min[i] = opt.x[i]
#
c[i] = -1
opt = scipy.optimize.linprog(c,A_ub=A,b_ub=b,bounds=(-inf,inf))
if not opt.success:
raise ValueError, errmsg.format(opt)
max[i] = opt.x[i]
return min,max
def aabb(self):
return Vector3(self._min), Vector3(self._max)
def center(self):
return .5*(self._min + self._max)
def dim(self):
return self._max - self._min
def __call__(self,pt,pad):
for p in self._inHalfSpaces:
if not p(pt,pad):
return False
return True
#####
## surface construction and manipulation
#####
def gtsSurface2Facets(surf,**kw):
"""Construct facets from given GTS surface. \*\*kw is passed to utils.facet."""
import gts
return [utils.facet([v.coords() for v in face.vertices()],**kw) for face in surf.faces()]
def sweptPolylines2gtsSurface(pts,threshold=0,capStart=False,capEnd=False):
"""Create swept suface (as GTS triangulation) given same-length sequences of points (as 3-tuples).
If threshold is given (>0), then
* degenerate faces (with edges shorter than threshold) will not be created
* gts.Surface().cleanup(threshold) will be called before returning, which merges vertices mutually closer than threshold. In case your pts are closed (last point concident with the first one) this will the surface strip of triangles. If you additionally have capStart==True and capEnd==True, the surface will be closed.
.. note:: capStart and capEnd make the most naive polygon triangulation (diagonals) and will perhaps fail for non-convex sections.
.. warning:: the algorithm connects points sequentially; if two polylines are mutually rotated or have inverse sense, the algorithm will not detect it and connect them regardless in their given order.
"""
import gts # will raise an exception in gts-less builds
if not len(set([len(pts1) for pts1 in pts]))==1: raise RuntimeError("Polylines must be all of the same length!")
vtxs=[[gts.Vertex(x,y,z) for x,y,z in pts1] for pts1 in pts]
sectEdges=[[gts.Edge(vtx[i],vtx[i+1]) for i in xrange(0,len(vtx)-1)] for vtx in vtxs]
interSectEdges=[[] for i in range(0,len(vtxs)-1)]
for i in range(0,len(vtxs)-1):
for j in range(0,len(vtxs[i])):
interSectEdges[i].append(gts.Edge(vtxs[i][j],vtxs[i+1][j]))
if j<len(vtxs[i])-1: interSectEdges[i].append(gts.Edge(vtxs[i][j],vtxs[i+1][j+1]))
if threshold>0: # replace edges of zero length with None; their faces will be skipped
def fixEdges(edges):
for i,e in enumerate(edges):
if (Vector3(e.v1.x,e.v1.y,e.v1.z)-Vector3(e.v2.x,e.v2.y,e.v2.z)).norm()<threshold: edges[i]=None
for ee in sectEdges: fixEdges(ee)
for ee in interSectEdges: fixEdges(ee)
surf=gts.Surface()
for i in range(0,len(vtxs)-1):
for j in range(0,len(vtxs[i])-1):
ee1=interSectEdges[i][2*j+1],sectEdges[i+1][j],interSectEdges[i][2*j]
ee2=sectEdges[i][j],interSectEdges[i][2*j+2],interSectEdges[i][2*j+1]
if None not in ee1: surf.add(gts.Face(interSectEdges[i][2*j+1],sectEdges[i+1][j],interSectEdges[i][2*j]))
if None not in ee2: surf.add(gts.Face(sectEdges[i][j],interSectEdges[i][2*j+2],interSectEdges[i][2*j+1]))
def doCap(vtx,edg,start):
ret=[]
eFan=[edg[0]]+[gts.Edge(vtx[i],vtx[0]) for i in range(2,len(vtx))]
for i in range(1,len(edg)):
ret+=[gts.Face(eFan[i-1],eFan[i],edg[i]) if start else gts.Face(eFan[i-1],edg[i],eFan[i])]
return ret
caps=[]
if capStart: caps+=doCap(vtxs[0],sectEdges[0],start=True)
if capEnd: caps+=doCap(vtxs[-1],sectEdges[-1],start=False)
for cap in caps: surf.add(cap)
if threshold>0: surf.cleanup(threshold)
return surf
def gtsSurfaceBestFitOBB(surf):
"""Return (Vector3 center, Vector3 halfSize, Quaternion orientation) describing
best-fit oriented bounding box (OBB) for the given surface. See cloudBestFitOBB
for details."""
import gts
pts=[Vector3(v.x,v.y,v.z) for v in surf.vertices()]
return cloudBestFitOBB(tuple(pts))
def revolutionSurfaceMeridians(sects,angles,origin=Vector3().Zero,orientation=Quaternion().Identity):
"""Revolution surface given sequences of 2d points and sequence of corresponding angles,
returning sequences of 3d points representing meridian sections of the revolution surface.
The 2d sections are turned around z-axis, but they can be transformed
using the origin and orientation arguments to give arbitrary orientation."""
import math
def toGlobal(x,y,z):
return tuple(origin+orientation*(Vector3(x,y,z)))
return [[toGlobal(x2d*math.cos(angles[i]),x2d*math.sin(angles[i]),y2d) for x2d,y2d in sects[i]] for i in range(0,len(sects))]
########
## packing generators
########
def regularOrtho(predicate,radius,gap,**kw):
"""Return set of spheres in regular orthogonal grid, clipped inside solid given by predicate.
Created spheres will have given radius and will be separated by gap space."""
ret=[]
mn,mx=predicate.aabb()
if(max([mx[i]-mn[i] for i in 0,1,2])==float('inf')): raise ValueError("Aabb of the predicate must not be infinite (didn't you use union | instead of intersection & for unbounded predicate such as notInNotch?");
xx,yy,zz=[arange(mn[i]+radius,mx[i]-radius,2*radius+gap) for i in 0,1,2]
for xyz in itertools.product(xx,yy,zz):
if predicate(xyz,radius): ret+=[utils.sphere(xyz,radius=radius,**kw)]
if (len(ret)==0):
warnings.warn('No spheres are produced by regularOrtho-function',category=RuntimeWarning)
return ret
def regularHexa(predicate,radius,gap,**kw):
"""Return set of spheres in regular hexagonal grid, clipped inside solid given by predicate.
Created spheres will have given radius and will be separated by gap space."""
ret=[]
a=2*radius+gap
hy,hz=a*sqrt(3)/2.,a*sqrt(6)/3.
mn,mx=predicate.aabb()
dim=[mx[i]-mn[i] for i in 0,1,2]
if(max(dim)==float('inf')): raise ValueError("Aabb of the predicate must not be infinite (didn't you use union | instead of intersection & for unbounded predicate such as notInNotch?");
ii,jj,kk=[range(0,int(dim[0]/a)+1),range(0,int(dim[1]/hy)+1),range(0,int(dim[2]/hz)+1)]
for i,j,k in itertools.product(ii,jj,kk):
#Simple HCP-lattice packing
#http://en.wikipedia.org/wiki/Close-packing_of_equal_spheres#Simple_hcp_lattice
coordSph = Vector3((2*i + ((j + k) % 2 ) ),
(sqrt(3.)*(j + 1./3.*(k % 2))),
(2.*sqrt(6.)/3.*k))*(a/2.0) + mn
if predicate(coordSph,radius): ret+=[utils.sphere(coordSph,radius=radius,**kw)]
if (len(ret)==0):
warnings.warn('No spheres are produced by regularHexa-function',category=RuntimeWarning)
return ret
def filterSpherePack(predicate,spherePack,returnSpherePack=None,**kw):
"""Using given SpherePack instance, return spheres that satisfy predicate.
It returns either a :yref:`yade._packSpheres.SpherePack` (if returnSpherePack) or a list.
The packing will be recentered to match the predicate and warning is given if the predicate
is larger than the packing."""
if returnSpherePack==None:
warnings.warn('The default behavior will change; specify returnSpherePack=True for the new behavior, and False to get rid of this warning (your code will break in the future, however). The returned SpherePack object can be added to the simulation using SpherePack.toSimulation()',category=FutureWarning)
returnSpherePack=False
mn,mx=predicate.aabb()
dimP,centP=predicate.dim(),predicate.center()
dimS,centS=spherePack.dim(),spherePack.center()
if dimP[0]>dimS[0] or dimP[1]>dimS[1] or dimP[2]>dimS[2]: warnings.warn("Packing's dimension (%s) doesn't fully contain dimension of the predicate (%s)."%(dimS,dimP))
spherePack.translate(centP-centS)
if returnSpherePack:
ret=SpherePack()
for c,r in spherePack:
if predicate(c,r): ret.add(c,r)
return ret
else:
# return particles to be added to O.bodies
ret=[]
for s in spherePack:
if predicate(s[0],s[1]): ret+=[utils.sphere(s[0],radius=s[1],**kw)]
return ret
def _memoizePacking(memoizeDb,sp,radius,rRelFuzz,wantPeri,fullDim,noPrint=False):
if not memoizeDb: return
import cPickle,sqlite3,time,os
if os.path.exists(memoizeDb):
conn=sqlite3.connect(memoizeDb)
else:
conn=sqlite3.connect(memoizeDb)
c=conn.cursor()
c.execute('create table packings (radius real, rRelFuzz real, dimx real, dimy real, dimz real, N integer, timestamp real, periodic integer, pack blob)')
c=conn.cursor()
packBlob=buffer(cPickle.dumps(sp.toList(),cPickle.HIGHEST_PROTOCOL))
packDim=sp.cellSize if wantPeri else fullDim
c.execute('insert into packings values (?,?,?,?,?,?,?,?,?)',(radius,rRelFuzz,packDim[0],packDim[1],packDim[2],len(sp),time.time(),wantPeri,packBlob,))
c.close()
conn.commit()
if not noPrint: print "Packing saved to the database",memoizeDb
def _getMemoizedPacking(memoizeDb,radius,rRelFuzz,x1,y1,z1,fullDim,wantPeri,fillPeriodic,spheresInCell,memoDbg=False,noPrint=False):
"""Return suitable SpherePack read from *memoizeDb* if found, None otherwise.
:param fillPeriodic: whether to fill fullDim by repeating periodic packing
:param wantPeri: only consider periodic packings
"""
import os,os.path,sqlite3,time,cPickle,sys
if memoDbg and not noPrint:
def memoDbgMsg(s): print s
else:
def memoDbgMsg(s): pass
if not memoizeDb or not os.path.exists(memoizeDb):
if memoizeDb: memoDbgMsg("Database %s does not exist."%memoizeDb)
return None
# find suitable packing and return it directly
conn=sqlite3.connect(memoizeDb); c=conn.cursor();
try:
c.execute('select radius,rRelFuzz,dimx,dimy,dimz,N,timestamp,periodic from packings order by N')
except sqlite3.OperationalError:
raise RuntimeError("ERROR: database `"+memoizeDb+"' not compatible with randomDensePack (corrupt, deprecated format or not a db created by randomDensePack)")
for row in c:
R,rDev,X,Y,Z,NN,timestamp,isPeri=row[0:8]; scale=radius/R
rDev*=scale; X*=scale; Y*=scale; Z*=scale
memoDbgMsg("Considering packing (radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created %s"%(R,.5*rDev,NN,X,Y,Z,"periodic" if isPeri else "non-periodic",scale,time.asctime(time.gmtime(timestamp))))
if not isPeri and wantPeri: memoDbgMsg("REJECT: is not periodic, which is requested."); continue
if wantPeri and (X/x1>0.9 or X/x1<0.6): memoDbgMsg("REJECT: initSize differs too much from scaled packing size."); continue
if (rRelFuzz==0 and rDev!=0) or (rRelFuzz!=0 and rDev==0) or (rRelFuzz!=0 and abs((rDev-rRelFuzz)/rRelFuzz)>1e-2): memoDbgMsg("REJECT: radius fuzz differs too much (%g, %g desired)"%(rDev,rRelFuzz)); continue # radius fuzz differs too much
if isPeri and wantPeri:
if spheresInCell>NN and spheresInCell>0: memoDbgMsg("REJECT: Number of spheres in the packing too small"); continue
if abs((y1/x1)/(Y/X)-1)>0.3 or abs((z1/x1)/(Z/X)-1)>0.3: memoDbgMsg("REJECT: proportions (y/x=%g, z/x=%g) differ too much from what is desired (%g, %g)."%(Y/X,Z/X,y1/x1,z1/x1)); continue
else:
if (X<fullDim[0] or Y<fullDim[1] or Z<fullDim[2]): memoDbgMsg("REJECT: not large enough"); continue # not large enough
memoDbgMsg("ACCEPTED");
if not noPrint: print "Found suitable packing in %s (radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created %s"%(memoizeDb,R,rDev,NN,X,Y,Z,"periodic" if isPeri else "non-periodic",scale,time.asctime(time.gmtime(timestamp)))
c.execute('select pack from packings where timestamp=?',(timestamp,))
sp=SpherePack(cPickle.loads(str(c.fetchone()[0])))
sp.scale(scale);
if isPeri and wantPeri:
sp.isPeriodic = True
sp.cellSize=(X,Y,Z);
if fillPeriodic: sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2]));
#sp.cellSize=(0,0,0) # resetting cellSize avoids warning when rotating
return sp
#if orientation: sp.rotate(*orientation.toAxisAngle())
#return filterSpherePack(predicate,sp,material=material)
#print "No suitable packing in database found, running",'PERIODIC compression' if wantPeri else 'triaxial'
#sys.stdout.flush()
def randomDensePack(predicate,radius,material=-1,dim=None,cropLayers=0,rRelFuzz=0.,spheresInCell=0,memoizeDb=None,useOBB=False,memoDbg=False,color=None,returnSpherePack=None):
"""Generator of random dense packing with given geometry properties, using TriaxialTest (aperiodic)
or PeriIsoCompressor (periodic). The periodicity depens on whether the spheresInCell parameter is given.
*O.switchScene()* magic is used to have clean simulation for TriaxialTest without deleting the original simulation.
This function therefore should never run in parallel with some code accessing your simulation.
:param predicate: solid-defining predicate for which we generate packing
:param spheresInCell: if given, the packing will be periodic, with given number of spheres in the periodic cell.
:param radius: mean radius of spheres
:param rRelFuzz: relative fuzz of the radius -- e.g. radius=10, rRelFuzz=.2, then spheres will have radii 10 ± (10*.2)), with an uniform distribution.
0 by default, meaning all spheres will have exactly the same radius.
:param cropLayers: (aperiodic only) how many layers of spheres will be added to the computed dimension of the box so that there no
(or not so much, at least) boundary effects at the boundaries of the predicate.
:param dim: dimension of the packing, to override dimensions of the predicate (if it is infinite, for instance)
:param memoizeDb: name of sqlite database (existent or nonexistent) to find an already generated packing or to store
the packing that will be generated, if not found (the technique of caching results of expensive computations
is known as memoization). Fuzzy matching is used to select suitable candidate -- packing will be scaled, rRelFuzz
and dimensions compared. Packing that are too small are dictarded. From the remaining candidate, the one with the
least number spheres will be loaded and returned.
:param useOBB: effective only if a inGtsSurface predicate is given. If true (not default), oriented bounding box will be
computed first; it can reduce substantially number of spheres for the triaxial compression (like 10× depending on
how much asymmetric the body is), see examples/gts-horse/gts-random-pack-obb.py
:param memoDbg: show packings that are considered and reasons why they are rejected/accepted
:param returnSpherePack: see the corresponding argument in :yref:`yade.pack.filterSpherePack`
:return: SpherePack object with spheres, filtered by the predicate.
"""
import sqlite3, os.path, cPickle, time, sys, _packPredicates, numpy
from math import pi
wantPeri=(spheresInCell>0)
if 'inGtsSurface' in dir(_packPredicates) and type(predicate)==inGtsSurface and useOBB:
center,dim,orientation=gtsSurfaceBestFitOBB(predicate.surf)
print "Best-fit oriented-bounding-box computed for GTS surface, orientation is",orientation
dim*=2 # gtsSurfaceBestFitOBB returns halfSize
else:
if not dim: dim=predicate.dim()
if max(dim)==float('inf'): raise RuntimeError("Infinite predicate and no dimension of packing requested.")
center=predicate.center()
orientation=None
if not wantPeri: fullDim=tuple([dim[i]+4*cropLayers*radius for i in 0,1,2])
else:
# compute cell dimensions now, as they will be compared to ones stored in the db
# they have to be adjusted to not make the cell to small WRT particle radius
fullDim=dim
cloudPorosity=0.25 # assume this number for the initial cloud (can be underestimated)
beta,gamma=fullDim[1]/fullDim[0],fullDim[2]/fullDim[0] # ratios β=y₀/x₀, γ=z₀/x₀
N100=spheresInCell/cloudPorosity # number of spheres for cell being filled by spheres without porosity
x1=radius*(1/(beta*gamma)*N100*(4/3.)*pi)**(1/3.)
y1,z1=beta*x1,gamma*x1; vol0=x1*y1*z1
maxR=radius*(1+rRelFuzz)
x1=max(x1,8*maxR); y1=max(y1,8*maxR); z1=max(z1,8*maxR); vol1=x1*y1*z1
N100*=vol1/vol0 # volume might have been increased, increase number of spheres to keep porosity the same
sp=_getMemoizedPacking(memoizeDb,radius,rRelFuzz,x1,y1,z1,fullDim,wantPeri,fillPeriodic=True,spheresInCell=spheresInCell,memoDbg=False)
if sp:
if orientation:
sp.cellSize=(0,0,0) # resetting cellSize avoids warning when rotating
sp.rotate(*orientation.toAxisAngle())
return filterSpherePack(predicate,sp,material=material,returnSpherePack=returnSpherePack)
else: print "No suitable packing in database found, running",'PERIODIC compression' if wantPeri else 'triaxial'
sys.stdout.flush()
O.switchScene(); O.resetThisScene() ### !!
if wantPeri:
# x1,y1,z1 already computed above
sp=SpherePack()
O.periodic=True
#O.cell.refSize=(x1,y1,z1)
O.cell.setBox((x1,y1,z1))
#print cloudPorosity,beta,gamma,N100,x1,y1,z1,O.cell.refSize
#print x1,y1,z1,radius,rRelFuzz
O.materials.append(FrictMat(young=3e10,density=2400))
num=sp.makeCloud(Vector3().Zero,O.cell.refSize,radius,rRelFuzz,spheresInCell,True)
O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=.05*radius),InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=5,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()
#print 'Resulting cellSize',sp.cellSize,'proportions',sp.cellSize[1]/sp.cellSize[0],sp.cellSize[2]/sp.cellSize[0]
# repetition to the required cell size will be done below, after memoizing the result
else:
assumedFinalDensity=0.6
V=(4.0/3.0)*pi*radius**3.0; N=assumedFinalDensity*fullDim[0]*fullDim[1]*fullDim[2]/V;
TriaxialTest(
numberOfGrains=int(N),radiusMean=radius,radiusStdDev=rRelFuzz,
# upperCorner is just size ratio, if radiusMean is specified
upperCorner=fullDim,
## no need to touch any the following
noFiles=True,lowerCorner=[0,0,0],sigmaIsoCompaction=-1e7,sigmaLateralConfinement=-1e5,compactionFrictionDeg=1,StabilityCriterion=.05,strainRate=.2,thickness=-1,maxWallVelocity=.1,wallOversizeFactor=1.5,autoUnload=True,autoCompressionActivation=False).load()
while ( numpy.isnan(utils.unbalancedForce()) or utils.unbalancedForce()>0.005 ) :
O.run(100,True)
sp=SpherePack(); sp.fromSimulation()
O.switchScene() ### !!
_memoizePacking(memoizeDb,sp,radius,rRelFuzz,wantPeri,fullDim)
if wantPeri: sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2]))
if orientation:
sp.cellSize=(0,0,0); # reset periodicity to avoid warning when rotating periodic packing
sp.rotate(*orientation.toAxisAngle())
return filterSpherePack(predicate,sp,material=material,color=color,returnSpherePack=returnSpherePack)
def randomPeriPack(radius,initSize,rRelFuzz=0.0,memoizeDb=None,noPrint=False):
"""Generate periodic dense packing.
A cell of initSize is stuffed with as many spheres as possible, then we run periodic compression with PeriIsoCompressor, just like with
randomDensePack.
:param radius: mean sphere radius
:param rRelFuzz: relative fuzz of sphere radius (equal distribution); see the same param for randomDensePack.
:param initSize: initial size of the periodic cell.
:return: SpherePack object, which also contains periodicity information.
"""
from math import pi
sp=_getMemoizedPacking(memoizeDb,radius,rRelFuzz,initSize[0],initSize[1],initSize[2],fullDim=Vector3(0,0,0),wantPeri=True,fillPeriodic=False,spheresInCell=-1,memoDbg=True,noPrint=noPrint)
if sp: return sp
O.switchScene(); O.resetThisScene()
sp=SpherePack()
O.periodic=True
#O.cell.refSize=initSize
O.cell.setBox(initSize)
sp.makeCloud(Vector3().Zero,O.cell.refSize,radius,rRelFuzz,-1,True)
O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=.05*radius),InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=20,keepProportions=True),NewtonIntegrator(damping=.8)]
O.materials.append(FrictMat(young=30e9,frictionAngle=.1,poisson=.3,density=1e3))
for s in sp: O.bodies.append(utils.sphere(s[0],s[1]))
O.dt=utils.PWaveTimeStep()
O.timingEnabled=True
O.run(); O.wait()
ret=SpherePack()
ret.fromSimulation()
_memoizePacking(memoizeDb,ret,radius,rRelFuzz,wantPeri=True,fullDim=Vector3(0,0,0),noPrint=noPrint) # fullDim unused
O.switchScene()
return ret
def hexaNet( radius, cornerCoord=[0,0,0], xLength=1., yLength=0.5, mos=0.08, a=0.04, b=0.04, startAtCorner=True, isSymmetric=False, **kw ):
"""Definition of the particles for a hexagonal wire net in the x-y-plane for the WireMatPM.
:param radius: radius of the particle
:param cornerCoord: coordinates of the lower left corner of the net
:param xLenght: net length in x-direction
:param yLenght: net length in y-direction
:param mos: mesh opening size (horizontal distance between the double twists)
:param a: length of double-twist
:param b: height of single wire section
:param startAtCorner: if true the generation starts with a double-twist at the lower left corner
:param isSymmetric: defines if the net is symmetric with respect to the y-axis
:return: set of spheres which defines the net (net) and exact dimensions of the net (lx,ly).
.. note:: This packing works for the WireMatPM only. The particles at the corner are always generated first. For examples on how to use this packing see examples/WireMatPM. In order to create the proper interactions for the net the interaction radius has to be adapted in the simulation.
"""
# check input dimension
if(xLength<mos): raise ValueError("xLength must be greater than mos!");
if(yLength<2*a+b): raise ValueError("yLength must be greater than 2*a+b!");
xstart = cornerCoord[0]
ystart = cornerCoord[1]
z = cornerCoord[2]
ab = a+b
# number of double twisted sections in y-direction and real length ly
ny = int( (yLength-a)/ab ) + 1
ly = ny*a+(ny-1)*b
jump=0
# number of sections in x-direction and real length lx
if isSymmetric:
nx = int( xLength/mos ) + 1
lx = (nx-1)*mos
if not startAtCorner:
nx+=-1
else:
nx = int( (xLength-0.5*mos)/mos ) + 1
lx = (nx-1)*mos+0.5*mos
net = []
# generate corner particles
if startAtCorner:
if (ny%2==0): # if ny even no symmetry in y-direction
net+=[utils.sphere((xstart,ystart+ly,z),radius=radius,**kw)] # upper left corner
if isSymmetric:
net+=[utils.sphere((xstart+lx,ystart+ly,z),radius=radius,**kw)] # upper right corner
else:
net+=[utils.sphere((xstart+lx,ystart,z),radius=radius,**kw)] # lower right corner
else: # if ny odd symmetry in y-direction
if not isSymmetric:
net+=[utils.sphere((xstart+lx,ystart,z),radius=radius,**kw)] # lower right corner
net+=[utils.sphere((xstart+lx,ystart+ly,z),radius=radius,**kw)] # upper right corner
jump=1
else: # do not start at corner
if (ny%2==0): # if ny even no symmetry in y-direction
net+=[utils.sphere((xstart,ystart,z),radius=radius,**kw)] # lower left corner
if isSymmetric:
net+=[utils.sphere((xstart+lx,ystart,z),radius=radius,**kw)] # lower right corner
else:
net+=[utils.sphere((xstart+lx,ystart+ly,z),radius=radius,**kw)] # upper right corner
else: # if ny odd symmetry in y-direction
net+=[utils.sphere((xstart,ystart,z),radius=radius,**kw)] # lower left corner
net+=[utils.sphere((xstart,ystart+ly,z),radius=radius,**kw)] # upper left corner
if isSymmetric:
net+=[utils.sphere((xstart+lx,ystart,z),radius=radius,**kw)] # lower right corner
net+=[utils.sphere((xstart+lx,ystart+ly,z),radius=radius,**kw)] # upper right corner
xstart+=0.5*mos
# generate other particles
if isSymmetric:
for i in range(ny):
y = ystart + i*ab
for j in range(nx):
x = xstart + j*mos
# add two particles of one vertical section (double-twist)
net+=[utils.sphere((x,y,z),radius=radius,**kw)]
net+=[utils.sphere((x,y+a,z),radius=radius,**kw)]
# set values for next section
xstart = xstart - 0.5*mos*pow(-1,i+jump)
nx = int(nx + 1*pow(-1,i+jump))
else:
for i in range(ny):
y = ystart + i*ab
for j in range(nx):
x = xstart + j*mos
# add two particles of one vertical section (double-twist)
net+=[utils.sphere((x,y,z),radius=radius,**kw)]
net+=[utils.sphere((x,y+a,z),radius=radius,**kw)]
# set values for next section
xstart = xstart - 0.5*mos*pow(-1,i+jump)
return [net,lx,ly]
|