yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01816
[Branch ~yade-dev/yade/trunk] Rev 1718: 1. PeriIsoCompressor: add keepProportions parameter (on by default), decrease global updates inte...
------------------------------------------------------------
revno: 1718
committer: Václav Šmilauer <vaclav@flux>
branch nick: trunk
timestamp: Sat 2009-08-22 17:05:20 +0200
message:
1. PeriIsoCompressor: add keepProportions parameter (on by default), decrease global updates interval.
2. Omega: startSimulationLoop warns if there is no simulationLoop (how can that happen? It somehow can)
3. SpherePack: resets periodicity (and warns) if being rotated; cellSize is writable now.
4. pack.triaxialPack: renamed to pack.randomDensePack; fix loads of stuff in that; rename radiusStDev parameter to rRelFuzz; add epydoc documentation; fix erroneous detection of inGtsSurface being defined
5. wrapper: lock renderer when adding list of bodies; Omega().resetThisWorld to reset current world only.
6. Fix scripts/test/{gts-random-pack.py,gts-random-obb.py}
renamed:
scripts/test/gts-triax-pack-obb.py => scripts/test/gts-random-pack-obb.py
scripts/test/gts-triax-pack.py => scripts/test/gts-random-pack.py
modified:
core/Omega.cpp
extra/PeriodicInsertionSortCollider.cpp
extra/PeriodicInsertionSortCollider.hpp
pkg/dem/DataClass/SpherePack.cpp
pkg/dem/DataClass/SpherePack.hpp
py/_packSpheres.cpp
py/pack.py
py/yadeWrapper/yadeWrapper.cpp
scripts/test/periodic-compress.py
scripts/test/gts-random-pack-obb.py
scripts/test/gts-random-pack.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.
=== modified file 'core/Omega.cpp'
--- core/Omega.cpp 2009-08-19 10:48:06 +0000
+++ core/Omega.cpp 2009-08-22 15:05:20 +0000
@@ -100,6 +100,7 @@
resetRootBody();
rootBodyAnother=shared_ptr<MetaBody>(new MetaBody);
timeInit();
+ createSimulationLoop();
}
void Omega::timeInit(){
@@ -124,6 +125,7 @@
void Omega::startSimulationLoop(){
+ if(!simulationLoop){ LOG_ERROR("No Omega::simulationLoop? Creating one (please report bug)."); createSimulationLoop(); }
if (simulationLoop && !simulationLoop->looping()){
simulationPauseDuration += microsec_clock::local_time()-msStartingPauseTime;
simulationLoop->start();
=== modified file 'extra/PeriodicInsertionSortCollider.cpp'
--- extra/PeriodicInsertionSortCollider.cpp 2009-08-12 08:27:41 +0000
+++ extra/PeriodicInsertionSortCollider.cpp 2009-08-22 15:05:20 +0000
@@ -74,7 +74,7 @@
}
#endif
if((pmn1!=pmx1) || (pmn2!=pmx2)){
- LOG_FATAL("Body #"<<(pmn1!=pmx1?id1:id2)<<" spans over half of the cell size!");
+ LOG_FATAL("Body #"<<(pmn1!=pmx1?id1:id2)<<" spans over half of the cell size "<<dim<<" (axis="<<axis<<", min="<<(pmn1!=pmx1?mn1:mn2)<<", max="<<(pmn1!=mn1?mx1:mx2)<<")");
throw runtime_error(__FILE__ ": Body larger than half of the cell size encountered.");
}
periods[axis]=(int)(pmn1-pmn2);
@@ -372,11 +372,14 @@
if(!b->boundingVolume) continue;
for(int i=0; i<3; i++) maxSpan=max(maxSpan,b->boundingVolume->max[i]-b->boundingVolume->min[i]);
}
+
}
if(maxDisplPerStep<0) maxDisplPerStep=1e-2*charLen; // this should be tuned somehow…
const long& step=rb->currentIteration;
Vector3r cellSize=rb->cellMax-rb->cellMin; //unused: Real cellVolume=cellSize[0]*cellSize[1]*cellSize[2];
Vector3r cellArea=Vector3r(cellSize[1]*cellSize[2],cellSize[0]*cellSize[2],cellSize[0]*cellSize[1]);
+ Real minSize=min(cellSize[0],min(cellSize[1],cellSize[2]));
+ if(minSize<2.1*maxSpan){ throw runtime_error("Minimum cell size is smaller than 2.1*span_of_the_biggest_body! (periodic collider requirement)"); }
if(((step%globalUpdateInt)==0) || avgStiffness<0 || sigma[0]<0 || sigma[1]<0 || sigma[2]<0){
Vector3r sumForces=Shop::totalForceInVolume(avgStiffness,rb);
sigma=Vector3r(sumForces[0]/cellArea[0],sumForces[1]/cellArea[1],sumForces[2]/cellArea[2]);
@@ -387,18 +390,28 @@
Vector3r cellGrow(Vector3r::ZERO);
// is the stress condition satisfied in all directions?
bool allStressesOK=true;
- Vector3r cg_;
- for(int axis=0; axis<3; axis++){
- // Δσ=ΔεE=(Δl/l)×(l×K/A) ↔ Δl=Δσ×A/K
- // FIXME: either NormalShearInteraction::{kn,ks} is computed wrong or we have dimensionality problem here
- // FIXME: that is why the fixup 1e-4 is needed here
- // FIXME: or perhaps maxDisplaPerStep=1e-2*charLen is too big??
- cellGrow[axis]=1e-4*(sigma[axis]-sigmaGoal)*cellArea[axis]/(avgStiffness>0?avgStiffness:1);
- if(abs(cellGrow[axis])>maxDisplPerStep) cellGrow[axis]=Mathr::Sign(cellGrow[axis])*maxDisplPerStep;
- cellGrow[axis]=max(cellGrow[axis],-(cellSize[0]-2.1*maxSpan));
- // crude way of predicting sigma, for steps when it is not computed from intrs
- if(avgStiffness>0) sigma[axis]-=cellGrow[axis]*avgStiffness;
- if(abs((sigma[axis]-sigmaGoal)/sigmaGoal)>5e-3) allStressesOK=false;
+ if(keepProportions){ // the same algo as below, but operating on quantitites averaged over all dimensions
+ Real sigAvg=(sigma[0]+sigma[1]+sigma[2])/3., avgArea=(cellArea[0]+cellArea[1]+cellArea[2])/3.;
+ Real grow=1e-4*(sigAvg-sigmaGoal)*avgArea/(avgStiffness>0?avgStiffness:1);
+ if(abs(grow)>maxDisplPerStep) grow=Mathr::Sign(grow)*maxDisplPerStep;
+ grow=max(grow,-(minSize-2.1*maxSpan));
+ if(avgStiffness>0) { sigma-=(grow*avgStiffness)*Vector3r::ONE; sigAvg-=grow*avgStiffness; }
+ if(abs((sigAvg-sigmaGoal)/sigmaGoal)>5e-3) allStressesOK=false;
+ cellGrow=Vector3r(grow,grow,grow);
+ }
+ else{ // handle each dimension separately
+ for(int axis=0; axis<3; axis++){
+ // Δσ=ΔεE=(Δl/l)×(l×K/A) ↔ Δl=Δσ×A/K
+ // FIXME: either NormalShearInteraction::{kn,ks} is computed wrong or we have dimensionality problem here
+ // FIXME: that is why the fixup 1e-4 is needed here
+ // FIXME: or perhaps maxDisplaPerStep=1e-2*charLen is too big??
+ cellGrow[axis]=1e-4*(sigma[axis]-sigmaGoal)*cellArea[axis]/(avgStiffness>0?avgStiffness:1);
+ if(abs(cellGrow[axis])>maxDisplPerStep) cellGrow[axis]=Mathr::Sign(cellGrow[axis])*maxDisplPerStep;
+ cellGrow[axis]=max(cellGrow[axis],-(cellSize[axis]-2.1*maxSpan));
+ // crude way of predicting sigma, for steps when it is not computed from intrs
+ if(avgStiffness>0) sigma[axis]-=cellGrow[axis]*avgStiffness;
+ if(abs((sigma[axis]-sigmaGoal)/sigmaGoal)>5e-3) allStressesOK=false;
+ }
}
TRVAR4(cellGrow,sigma,sigmaGoal,avgStiffness);
rb->cellMin-=.5*cellGrow; rb->cellMax+=.5*cellGrow;
=== modified file 'extra/PeriodicInsertionSortCollider.hpp'
--- extra/PeriodicInsertionSortCollider.hpp 2009-08-12 08:27:41 +0000
+++ extra/PeriodicInsertionSortCollider.hpp 2009-08-22 15:05:20 +0000
@@ -166,9 +166,11 @@
size_t state;
//! python command to be run when reaching the last specified stress
string doneHook;
+ //! whether to exactly keep proportions of the cell
+ bool keepProportions;
void action(MetaBody*);
- PeriIsoCompressor(): avgStiffness(-1), maxDisplPerStep(-1), sumForces(Vector3r::ZERO), sigma(Vector3r::ZERO), charLen(-1), maxSpan(-1), maxUnbalanced(1e-4), currUnbalanced(-1), globalUpdateInt(100), state(0){}
- REGISTER_ATTRIBUTES(StandAloneEngine,(stresses)(charLen)(maxUnbalanced)(globalUpdateInt)(state)(doneHook)(maxSpan));
+ PeriIsoCompressor(): avgStiffness(-1), maxDisplPerStep(-1), sumForces(Vector3r::ZERO), sigma(Vector3r::ZERO), charLen(-1), maxSpan(-1), maxUnbalanced(1e-4), currUnbalanced(-1), globalUpdateInt(20), state(0), keepProportions(true){}
+ REGISTER_ATTRIBUTES(StandAloneEngine,(stresses)(charLen)(maxUnbalanced)(globalUpdateInt)(state)(doneHook)(maxSpan)(keepProportions));
DECLARE_LOGGER;
REGISTER_CLASS_AND_BASE(PeriIsoCompressor,StandAloneEngine);
};
=== modified file 'pkg/dem/DataClass/SpherePack.cpp'
--- pkg/dem/DataClass/SpherePack.cpp 2009-08-12 08:27:41 +0000
+++ pkg/dem/DataClass/SpherePack.cpp 2009-08-22 15:05:20 +0000
@@ -115,6 +115,7 @@
void SpherePack::cellFill(Vector3r vol){
Vector3<int> count;
for(int i=0; i<3; i++) count[i]=(int)(ceil(vol[i]/cellSize[i]));
+ LOG_DEBUG("Filling volume "<<vol<<" with cell "<<cellSize<<", repeat counts are "<<count);
cellRepeat(count);
}
=== modified file 'pkg/dem/DataClass/SpherePack.hpp'
--- pkg/dem/DataClass/SpherePack.hpp 2009-08-12 20:12:21 +0000
+++ pkg/dem/DataClass/SpherePack.hpp 2009-08-22 15:05:20 +0000
@@ -83,8 +83,11 @@
// transformations
void translate(const Vector3r& shift){ FOREACH(Sph& s, pack) s.c+=shift; }
- void rotate(const Vector3r& axis, Real angle){ Vector3r mid=midPt(); Quaternionr q(axis,angle); FOREACH(Sph& s, pack) s.c=q*(s.c-mid)+mid; }
- void scale(Real scale){ Vector3r mid=midPt(); FOREACH(Sph& s, pack) {s.c=scale*(s.c-mid)+mid; s.r*=abs(scale); } }
+ void rotate(const Vector3r& axis, Real angle){
+ if(cellSize!=Vector3r::ZERO) { LOG_WARN("Periodicity reset when rotating periodic packing."); cellSize=Vector3r::ZERO; }
+ Vector3r mid=midPt(); Quaternionr q(axis,angle); FOREACH(Sph& s, pack) s.c=q*(s.c-mid)+mid;
+ }
+ void scale(Real scale){ Vector3r mid=midPt(); cellSize*=scale; FOREACH(Sph& s, pack) {s.c=scale*(s.c-mid)+mid; s.r*=abs(scale); } }
// iteration
#ifdef YADE_PYTHON
=== modified file 'py/_packSpheres.cpp'
--- py/_packSpheres.cpp 2009-08-11 10:37:21 +0000
+++ py/_packSpheres.cpp 2009-08-22 15:05:20 +0000
@@ -15,7 +15,7 @@
.def("aabb",&SpherePack::aabb_py,"Get axis-aligned bounding box coordinates, as 2 3-tuples.")
.def("dim",&SpherePack::dim,"Return dimensions of the packing in terms of aabb(), as a 3-tuple.")
.def("center",&SpherePack::midPt,"Return coordinates of the bounding box center.")
- .def_readonly("cellSize",&SpherePack::cellSize,"Size of periodic cell; is Vector3(0,0,0) if not periodic.")
+ .def_readwrite("cellSize",&SpherePack::cellSize,"Size of periodic cell; is Vector3(0,0,0) if not periodic. (Change this property only if you know what you're doing).")
.def("cellFill",&SpherePack::cellFill,"Repeat the packing (if periodic) so that the results has dim() >= given size. The packing retains periodicity, but changes cellSize. Raises exception for non-periodic packing.")
.def("cellRepeat",&SpherePack::cellRepeat,"Repeat the packing given number of times in each dimension. Periodicity is retained, cellSize changes. Raises exception for non-periodic packing.")
.def("relDensity",&SpherePack::relDensity,"Relative packing density, measured as sum of spheres' volumes / aabb volume.\n(Sphere overlaps are ignored.)")
=== modified file 'py/pack.py'
--- py/pack.py 2009-08-21 14:14:26 +0000
+++ py/pack.py 2009-08-22 15:05:20 +0000
@@ -24,7 +24,7 @@
except ImportError: pass
# make c++ predicates available in this module
-from _packPredicates import *
+from _packPredicates import * ## imported in randomDensePack as well
# import SpherePack
from _packSpheres import *
from _packObb import *
@@ -202,31 +202,40 @@
if predicate(s[0],s[1]): ret+=[utils.sphere(s[0],radius=s[1],**kw)]
return ret
-def triaxialPack(predicate,radius,dim=None,cropLayers=0,radiusStDev=0.,assumedFinalDensity=.6,spheresInCell=0,memoizeDb=None,useOBB=True,**kw):
- """Generator of triaxial packing, using TriaxialTest. Radius is radius of spheres, radiusStDev is its standard deviation.
- By default, all spheres are of the same radius. cropLayers is 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.
- assumedFinalDensity should be OK as it is, it is used to compute necessary number of spheres for the packing.
-
- The memoizeDb parameter can be passed a file (existent or nonexistent). If the file exists, it will be first looked
- for a suitable packing that was previously saved already (known as memoization). Saved packings will be scaled to
- requested sphere radius; those that are smaller are distcarded as well as those with different radiusStDev. From
- the remaining ones, the one with the least spheres will be loaded and returned. If no suitable packing is found, it
- is generated as usually, but saved into the database for later use.
-
- useOBB is effective only if a inGtsSurface predicate is given. If true (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 scripts/test/gts-triax-pack-obb.py.
-
- O.switchWorld() magic is used to have clean simulation for TriaxialTest without deleting the original simulation.
+def randomDensePack(predicate,radius,dim=None,cropLayers=0,rRelFuzz=0.,spheresInCell=0,memoizeDb=None,useOBB=True,memoDbg=False,**kw):
+ """Generator of random dense packing with given geometry properties, using TriaxialTest (aperiodic)
+ or PeriIsoCompressor (periodic). The priodicity depens on whether the spheresInCell parameter is given.
+
+ L{O.switchWorld()<Omega.switchWorld>} 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)).
+ 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 (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 scripts/test/gts-triax-pack-obb.py.
+ @param memoDbg: show packigns that are considered and reasons why they are rejected/accepted
+
+ @return: SpherePack object with spheres, filtered by the predicate.
"""
- import sqlite3, os.path, cPickle, time, sys
+ import sqlite3, os.path, cPickle, time, sys, _packPredicates
from yade import log
from math import pi
wantPeri=(spheresInCell>0)
- if 'inGtsSurface' in dir() and type(predicate)==inGtsSurface and useOBB:
+ 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()
@@ -234,52 +243,68 @@
center=predicate.center()
orientation=None
if not wantPeri: fullDim=tuple([dim[i]+4*cropLayers*radius for i in 0,1,2])
- else: fullDim=dim
+ 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
+ maxR=radius*(1+rRelFuzz)
+ x1=max(x1,8*maxR); y1=max(y1,8*maxR); z1=max(z1,8*maxR) # this might make the packing looser, oh well...
if(memoizeDb and os.path.exists(memoizeDb)):
+ if memoDbg:
+ def memoDbgMsg(s): print s
+ else:
+ def memoDbgMsg(s): pass
# find suitable packing and return it directly
conn=sqlite3.connect(memoizeDb); c=conn.cursor();
- c.execute('select radius,radiusStDev,dimx,dimy,dimz,N,timestamp,periodic from packings order by N')
+ try:
+ c.execute('select radius,rRelFuzz,dimx,dimy,dimz,N,timestamp,periodic from packings order by N')
+ except OperationalError:
+ raise RuntimeError("ERROR: database",memoizeDb," not compatible with randomDensePack (deprecated format or not 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
- print "Considering packing (radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created %s"%(R,rDev,NN,X,Y,Z,"periodic" if isPeri else "non-periodic",scale,time.asctime(time.gmtime(timestamp)))
- if (radiusStDev==0 and rDev!=0) or (radiusStDev==0 and rDev!=0) or (radiusStDev!=0 and abs((rDev-radiusStDev)/radiusStDev)>1e-2): continue # radius fuzz differs too much
+ 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 (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: continue
- if abs((fullDim[1]/fullDim[0])/(dimy/dimx)-1)>0.2 or abs((fullDim[2]/fullDim[0])/(dimz/dimx)-1)>0.2: continue # proportions differing too much
+ if spheresInCell>NN: memoDbgMsg("REJECT: Number of spheres in the packing too small"); continue
+ if abs((fullDim[1]/fullDim[0])/(Y/X)-1)>0.3 or abs((fullDim[2]/fullDim[0])/(Z/X)-1)>0.3: memoDbgMsg("REJECT: proportions differ too much from what is desired."); continue
else:
- if (X<fullDim[0] or Y<fullDim[1] or Z<fullDim[2]): continue # not large enough
- print "Found suitable packing in database (radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created %s"%(R,rDev,NN,X,Y,Z,"periodic" if isPeri else "non-periodic",scale,time.asctime(time.gmtime(timestamp)))
+ if (X<fullDim[0] or Y<fullDim[1] or Z<fullDim[2]): memoDbgMsg("REJECT: not large enough"); continue # not large enough
+ memoDbgMsg("ACCEPTED");
+ 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])))
+ if isPeri and wantPeri:
+ sp.cellSize=(X,Y,Z); sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2])); sp.cellSize=(0,0,0) # resetting cellSize avoids warning when rotating
sp.scale(scale);
if orientation: sp.rotate(*orientation.ToAxisAngle())
return filterSpherePack(predicate,sp,**kw)
- print "No suitable packing in database found, running triaxial"
+ print "No suitable packing in database found, running",'PERIODIC compression' if wantPeri else 'triaxial'
sys.stdout.flush()
- O.switchWorld() ### !!
+ O.switchWorld(); O.resetThisWorld() ### !!
if wantPeri:
- #O.reset() # doesn't (shouldn't) affect the original simulation
+ # x1,y1,z1 already computed above
sp=SpherePack()
- 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
O.periodicCell=((0,0,0),(x1,y1,z1))
- print cloudPorosity,beta,gamma,N100,x1,y1,z1,O.periodicCell
- num=sp.makeCloud(O.periodicCell[0],O.periodicCell[1],radius,radiusStDev,spheresInCell,True)
- O.engines=[BexResetter(),BoundingVolumeMetaEngine([InteractingSphere2AABB()]),PeriodicInsertionSortCollider(),InteractionDispatchers([ef2_Sphere_Sphere_Dem3DofGeom()],[SimpleElasticRelationships()],[Law2_Dem3Dof_Elastic_Elastic()]),PeriIsoCompressor(charLen=radius/5.,stresses=[100e9,1e9],maxUnbalanced=1e-2,doneHook='O.pause();'),NewtonsDampedLaw(damping=.6)]
+ #print cloudPorosity,beta,gamma,N100,x1,y1,z1,O.periodicCell
+ #print x1,y1,z1,radius,rRelFuzz
+ num=sp.makeCloud(O.periodicCell[0],O.periodicCell[1],radius,rRelFuzz,spheresInCell,True)
+ O.engines=[BexResetter(),BoundingVolumeMetaEngine([InteractingSphere2AABB()]),PeriodicInsertionSortCollider(),InteractionDispatchers([ef2_Sphere_Sphere_Dem3DofGeom()],[SimpleElasticRelationships()],[Law2_Dem3Dof_Elastic_Elastic()]),PeriIsoCompressor(charLen=radius/5.,stresses=[100e9,1e9],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=5),NewtonsDampedLaw(damping=.6)]
for s in sp: O.bodies.append(utils.sphere(s[0],s[1],density=1000))
O.dt=utils.PWaveTimeStep()
- #for i in range(10): O.step()
O.run(); O.wait()
sp=SpherePack(); sp.fromSimulation()
- sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2]))
+ # repetition to the required cell size will be done below, after memoizing the result
else:
+ assumedFinalDensity=0.6
V=(4/3)*pi*radius**3; N=assumedFinalDensity*fullDim[0]*fullDim[1]*fullDim[2]/V;
TriaxialTest(
- numberOfGrains=int(N),radiusMean=radius,radiusStdDev=radiusStDev,
+ 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
@@ -294,15 +319,21 @@
else:
conn=sqlite3.connect(memoizeDb)
c=conn.cursor()
- c.execute('create table packings (radius real, radiusStDev real, dimx real, dimy real, dimz real, N integer, timestamp real, periodic integer, pack blob)')
+ 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_pointsAsTuples(),cPickle.HIGHEST_PROTOCOL))
- c.execute('insert into packings values (?,?,?,?,?,?,?,?,?)',(radius,radiusStDev,fullDim[0],fullDim[1],fullDim[2],len(sp),time.time(),wantPeri,packBlob,))
+ 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()
print "Packing saved to the database",memoizeDb
+ if wantPeri: sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2]))
if orientation: sp.rotate(*orientation.ToAxisAngle())
return filterSpherePack(predicate,sp,**kw)
+# compatibility with the deprecated name, can be removed in the future
+def triaxialPack(*args,**kw):
+ import warnings; warnings.warn("pack.triaxialPack was renamed to pack.randomDensePack, update your code!",DeprecationWarning,stacklevel=2);
+ return randomDensePack(*args,**kw)
=== modified file 'py/yadeWrapper/yadeWrapper.cpp'
--- py/yadeWrapper/yadeWrapper.cpp 2009-08-21 11:44:00 +0000
+++ py/yadeWrapper/yadeWrapper.cpp 2009-08-22 15:05:20 +0000
@@ -134,6 +134,21 @@
}
body_id_t insert(shared_ptr<Body> b){ return proxee->insert(b); }
vector<body_id_t> insertList(vector<shared_ptr<Body> > bb){
+ /* prevent crash when adding lots of bodies (not clear why it happens exactly, bt is like this:
+
+ #3 <signal handler called>
+ #4 0x000000000052483f in boost::detail::atomic_increment (pw=0x8089) at /usr/include/boost/detail/sp_counted_base_gcc_x86.hpp:66
+ #5 0x00000000005248b3 in boost::detail::sp_counted_base::add_ref_copy (this=0x8081) at /usr/include/boost/detail/sp_counted_base_gcc_x86.hpp:133
+ #6 0x00000000005249ca in shared_count (this=0x7fff2e44db48, r=@0x7f08ffd692b8) at /usr/include/boost/detail/shared_count.hpp:227
+ #7 0x00000000005258e3 in shared_ptr (this=0x7fff2e44db40) at /usr/include/boost/shared_ptr.hpp:165
+ #8 0x0000000000505cff in BodyRedirectionVectorIterator::getValue (this=0x846f040) at /home/vaclav/yade/trunk/core/containers/BodyRedirectionVector.cpp:47
+ #9 0x00007f0908af41ce in BodyContainerIteratorPointer::operator* (this=0x7fff2e44db60) at /home/vaclav/yade/build-trunk/include/yade-trunk/yade/core/BodyContainer.hpp:63
+ #10 0x00007f0908af420a in boost::foreach_detail_::deref<BodyContainer, mpl_::bool_<false> > (cur=@0x7fff2e44db60) at /usr/include/boost/foreach.hpp:750
+ #11 0x00007f0908adc5a9 in OpenGLRenderingEngine::renderGeometricalModel (this=0x77f1240, rootBody=@0x1f49220) at pkg/common/RenderingEngine/OpenGLRenderingEngine/OpenGLRenderingEngine.cpp:441
+ #12 0x00007f0908adfb84 in OpenGLRenderingEngine::render (this=0x77f1240, rootBody=@0x1f49220, selection=-1) at pkg/common/RenderingEngine/OpenGLRenderingEngine/OpenGLRenderingEngine.cpp:232
+
+ */
+ boost::mutex::scoped_lock lock(Omega::instance().renderMutex);
vector<body_id_t> ret; FOREACH(shared_ptr<Body>& b, bb){ret.push_back(insert(b));} return ret;
}
python::tuple insertClump(vector<shared_ptr<Body> > bb){/*clump: first add constitutents, then add clump, then add constitutents to the clump, then update clump props*/
@@ -340,6 +355,7 @@
}
void reset(){Py_BEGIN_ALLOW_THREADS; OMEGA.reset(); Py_END_ALLOW_THREADS; }
+ void resetThisWorld(){Py_BEGIN_ALLOW_THREADS; OMEGA.joinSimulationLoop(); Py_END_ALLOW_THREADS; OMEGA.resetRootBody(); OMEGA.createSimulationLoop();}
void resetTime(){ OMEGA.getRootBody()->currentIteration=0; OMEGA.getRootBody()->simulationTime=0; OMEGA.timeInit(); }
void switchWorld(){ std::swap(OMEGA.rootBody,OMEGA.rootBodyAnother); }
@@ -631,7 +647,8 @@
.def("pause",&pyOmega::pause,"Stop simulation execution.\n(may be called from within the loop, and it will stop after the current step).")
.def("step",&pyOmega::step,"Advance the simulation by one step. Returns after the step will have finished.")
.def("wait",&pyOmega::wait,"Don't return until the simulation will have been paused. (Returns immediately if not running).")
- .def("reset",&pyOmega::reset,"Reset simulation completely.")
+ .def("reset",&pyOmega::reset,"Reset simulations completely (including another world!).")
+ .def("resetThisWorld",&pyOmega::resetThisWorld,"Reset current world.")
.def("switchWorld",&pyOmega::switchWorld,"Switch to alternative simulation (while keeping the old one). Calling the function again switches back to the first one. Note that most variables from the first simulation will still refer to the first simulation even after the switch\n(e.g. b=O.bodies[4]; O.switchWorld(); [b still refers to the body in the first simulation here])")
.def("labeledEngine",&pyOmega::labeled_engine_get,"Return instance of engine/functor with the given label. This function shouldn't be called by the user directly; every ehange in O.engines will assign respective global python variables according to labels.\n For example: \n\tO.engines=[InsertionSortCollider(label='collider')]\n\tcollider['nBins']=5 ## collider has become a variable after assignment to O.engines automatically)")
.def("resetTime",&pyOmega::resetTime,"Reset simulation time: step number, virtual and real time. (Doesn't touch anything else, including timings).")
=== renamed file 'scripts/test/gts-triax-pack-obb.py' => 'scripts/test/gts-random-pack-obb.py'
--- scripts/test/gts-triax-pack-obb.py 2009-07-17 20:50:55 +0000
+++ scripts/test/gts-random-pack-obb.py 2009-08-22 15:05:20 +0000
@@ -12,9 +12,9 @@
# fill this solid with triaxial packing; it will compute minimum-volume oriented bounding box
# to minimize the number of throw-away spheres.
# It does away with about 3k spheres for radius 3e-2
-O.bodies.append(pack.triaxialPack(pack.inGtsSurface(surf),radius=3e-2,radiusStDev=1e-1,memoizeDb='/tmp/gts-triax-packings.sqlite'))
+O.bodies.append(pack.randomDensePack(pack.inGtsSurface(surf),radius=3e-2,rRelFuzz=1e-1,memoizeDb='/tmp/gts-triax-packings.sqlite'))
# translate the surface away and pack it again with sphere, but without the oriented bounding box (useOBB=False)
# Here, we need 20k spheres (with more or less the same result)
surf.translate(0,0,1);
O.bodies.append(pack.gtsSurface2Facets(surf,color=(1,0,0)))
-O.bodies.append(pack.triaxialPack(pack.inGtsSurface(surf),radius=3e-2,radiusStDev=1e-1,memoizeDb='/tmp/gts-triax-packings.sqlite',useOBB=False))
+O.bodies.append(pack.randomDensePack(pack.inGtsSurface(surf),radius=3e-2,rRelFuzz=1e-1,memoizeDb='/tmp/gts-triax-packings.sqlite',useOBB=False))
=== renamed file 'scripts/test/gts-triax-pack.py' => 'scripts/test/gts-random-pack.py'
--- scripts/test/gts-triax-pack.py 2009-08-12 08:27:41 +0000
+++ scripts/test/gts-random-pack.py 2009-08-22 15:05:20 +0000
@@ -25,24 +25,22 @@
# 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
-# with given radius and standard deviation (see documentation of pack.triaxialPack)
+# with given radius and standard deviation (see documentation of pack.randomDensePack)
#
# 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!
memoizeDb='/tmp/gts-triax-packings.sqlite'
-O.bodies.append(pack.triaxialPack(pack.inGtsSurface(surf),radius=5e-3,radiusStDev=1e-4,memoizeDb=memoizeDb))
+O.bodies.append(pack.randomDensePack(pack.inGtsSurface(surf),radius=5e-3,rRelFuzz=1e-4,memoizeDb=memoizeDb))
# 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
if 1:
import gts
- horse=gts.read(open('horse.coarse.gts'))
- #horse.scale(.25,.25,.25)
- O.bodies.append(pack.gtsSurface2Facets(horse))
- O.bodies.append(pack.triaxialPack(pack.inGtsSurface(horse),radius=5e-3,memoizeDb=memoizeDb))
- ## NB: periodic triaxial doesn't work yet
- #horse.translate(0,.2,0)
- #O.bodies.append(pack.gtsSurface2Facets(horse))
- #O.bodies.append(pack.triaxialPack(pack.inGtsSurface(horse),radius=1,spheresInCell=2000))#,memoizeDb=memoizeDb))
+ horse=gts.read(open('horse.coarse.gts')) #; horse.scale(.25,.25,.25)
+ O.bodies.append(pack.gtsSurface2Facets(horse))
+ O.bodies.append(pack.randomDensePack(pack.inGtsSurface(horse),radius=5e-3,memoizeDb=memoizeDb))
+ horse.translate(.07,0,0)
+ O.bodies.append(pack.gtsSurface2Facets(horse))
+ O.bodies.append(pack.randomDensePack(pack.inGtsSurface(horse),radius=1e-3,spheresInCell=2000,memoizeDb=memoizeDb))
=== modified file 'scripts/test/periodic-compress.py'
--- scripts/test/periodic-compress.py 2009-08-12 08:27:41 +0000
+++ scripts/test/periodic-compress.py 2009-08-22 15:05:20 +0000
@@ -1,4 +1,4 @@
-O.periodicCell=(0,0,0),(20,20,20)
+O.periodicCell=(0,0,0),(20,20,10)
from yade import pack,log,timing
p=pack.SpherePack()
p.makeCloud(O.periodicCell[0],O.periodicCell[1],1,.5,700,True)
@@ -16,7 +16,7 @@
[SimpleElasticRelationships()],
[Law2_Dem3Dof_Elastic_Elastic()],
),
- PeriIsoCompressor(charLen=.5,stresses=[50e9,1e8],doneHook="print 'FINISHED'; O.pause() "),
+ PeriIsoCompressor(charLen=.5,stresses=[50e9,1e8],doneHook="print 'FINISHED'; O.pause() ",keepProportions=True),
NewtonsDampedLaw(damping=.4)
]
O.dt=utils.PWaveTimeStep()
