yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06122
[Branch ~yade-dev/yade/trunk] Rev 2537: 1. Implement an initial version of generating loose clump packings of predefined configuration wi...
------------------------------------------------------------
revno: 2537
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Tue 2010-11-09 09:27:14 +0100
message:
1. Implement an initial version of generating loose clump packings of predefined configuration with SpherePack (TODO: do not recompute clump properties when adding them using toSimulation for each of them separately, but use those precomputed on given configurations). See scripts/test/clumpPack.py
2. Make selection of clump member highlight the whole clump in OpenGL
2. Add Dong2010 article to articles done with yade (authors contacted for fulltext)
added:
scripts/test/clumpPack.py
modified:
doc/yade-articles.bib
pkg/common/OpenGLRenderer.cpp
pkg/common/OpenGLRenderer.hpp
pkg/dem/SpherePack.cpp
pkg/dem/SpherePack.hpp
py/pack/_packSpheres.cpp
py/pack/pack.py
py/wrapper/customConverters.cpp
--
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 'doc/yade-articles.bib'
--- doc/yade-articles.bib 2010-11-03 15:50:02 +0000
+++ doc/yade-articles.bib 2010-11-09 08:27:14 +0000
@@ -280,3 +280,18 @@
note = "in press"
}
+@article{Dang2010,
+ author = {H. K. Dang and M. A. Meguid},
+ title = {Algorithm to Generate a Discrete Element Specimen with Predefined Properties},
+ publisher = {ASCE},
+ year = {2010},
+ journal = {International Journal of Geomechanics},
+ volume = {10},
+ number = {2},
+ pages = {85-91},
+ keywords = {Discrete elements; Geotechnical engineering; Grain size; Porosity; Algorithms; Granular materials},
+ doi = "10.1061/(ASCE)GM.1943-5622.0000028"
+}
+
+
+
=== modified file 'pkg/common/OpenGLRenderer.cpp'
--- pkg/common/OpenGLRenderer.cpp 2010-11-02 12:02:13 +0000
+++ pkg/common/OpenGLRenderer.cpp 2010-11-09 08:27:14 +0000
@@ -134,7 +134,7 @@
if(!initDone) init();
assert(initDone);
- current_selection = selection;
+ selId = selection;
scene=_scene;
@@ -229,7 +229,7 @@
if(!b) continue;
if(b->shape && ((b->getGroupMask() & mask) || b->getGroupMask()==0)){
if(!id && b->state->blockedDOFs==0) continue;
- if(current_selection==b->getId()){glLightModelfv(GL_LIGHT_MODEL_AMBIENT,ambientColorSelected);}
+ if(selId==b->getId()){glLightModelfv(GL_LIGHT_MODEL_AMBIENT,ambientColorSelected);}
{ // write text
glColor3f(1.0-bgColor[0],1.0-bgColor[1],1.0-bgColor[2]);
unsigned d = b->state->blockedDOFs;
@@ -239,11 +239,11 @@
if(dof && id) sId += " ";
if(id) str += sId;
if(dof) str += sDof;
- const Vector3r& h(current_selection==b->getId() ? highlightEmission0 : Vector3r(1,1,1));
+ const Vector3r& h(selId==b->getId() ? highlightEmission0 : Vector3r(1,1,1));
glColor3v(h);
GLUtils::GLDrawText(str,bodyDisp[b->id].pos,h);
}
- if(current_selection == b->getId()){glLightModelfv(GL_LIGHT_MODEL_AMBIENT,ambientColorUnselected);}
+ if(selId == b->getId()){glLightModelfv(GL_LIGHT_MODEL_AMBIENT,ambientColorUnselected);}
}
}
}
@@ -323,14 +323,16 @@
// ignored in non-selection mode, use it always
glPushName(b->id);
+ bool highlight=(b->id==selId || b->clumpId==selId || b->shape->highlight);
glPushMatrix();
AngleAxisr aa(ori);
glTranslatef(pos[0],pos[1],pos[2]);
glRotatef(aa.angle()*Mathr::RAD_TO_DEG,aa.axis()[0],aa.axis()[1],aa.axis()[2]);
- if(current_selection==b->getId() || b->shape->highlight){
+ if(highlight){
// set hightlight
- const Vector3r& h(current_selection==b->getId() ? highlightEmission0 : highlightEmission1);
+ // different color for body highlighted by selection and by the shape attribute
+ const Vector3r& h((selId==b->id||selId==b->clumpId) ? highlightEmission0 : highlightEmission1);
glMaterialv(GL_FRONT_AND_BACK,GL_EMISSION,h);
glMaterialv(GL_FRONT_AND_BACK,GL_SPECULAR,h);
shapeDispatcher(b->shape,b->state,wire || b->shape->wire,viewInfo);
@@ -342,7 +344,7 @@
shapeDispatcher(b->shape,b->state,wire || b->shape->wire,viewInfo);
}
glPopMatrix();
- if(current_selection==b->getId() || b->shape->highlight){
+ if(highlight){
if(!b->bound || wire || b->shape->wire) GLUtils::GLDrawInt(b->getId(),pos);
else {
// move the label towards the camera by the bounding box so that it is not hidden inside the body
=== modified file 'pkg/common/OpenGLRenderer.hpp'
--- pkg/common/OpenGLRenderer.hpp 2010-10-13 16:23:08 +0000
+++ pkg/common/OpenGLRenderer.hpp 2010-11-09 08:27:14 +0000
@@ -23,8 +23,6 @@
public:
int _nothing; // remove later, only target for deprecated selectBodyLimit
- Body::id_t current_selection;
-
static const int numClipPlanes=3;
bool pointClipped(const Vector3r& p);
@@ -104,6 +102,7 @@
((bool,intrGeom,false,,"Render :yref:`Interaction::geom` objects."))
((bool,intrPhys,false,,"Render :yref:`Interaction::phys` objects"))
((int,mask,((void)"draw everything",~0),,"Bitmask for showing only bodies where ((mask & :yref:`Body::mask`)!=0)"))
+ ((Body::id_t,selId,Body::ID_NONE,,"Id of particle that was selected by the user."))
((vector<Se3r>,clipPlaneSe3,vector<Se3r>(numClipPlanes,Se3r(Vector3r::Zero(),Quaternionr::Identity())),,"Position and orientation of clipping planes"))
((vector<bool>,clipPlaneActive,vector<bool>(numClipPlanes,false),,"Activate/deactivate respective clipping planes"))
((vector<shared_ptr<GlExtraDrawer> >,extraDrawers,,,"Additional rendering components (:yref:`GlExtraDrawer`)."))
=== modified file 'pkg/dem/SpherePack.cpp'
--- pkg/dem/SpherePack.cpp 2010-10-30 17:08:51 +0000
+++ pkg/dem/SpherePack.cpp 2010-11-09 08:27:14 +0000
@@ -19,38 +19,36 @@
using namespace std;
using namespace boost;
-
-
-void SpherePack::fromList(const python::list& l){
+namespace py=boost::python;
+
+
+void SpherePack::fromList(const py::list& l){
pack.clear();
- size_t len=python::len(l);
+ size_t len=py::len(l);
for(size_t i=0; i<len; i++){
- const python::tuple& t=python::extract<python::tuple>(l[i]);
- python::extract<Vector3r> vec(t[0]);
- if(vec.check()) { pack.push_back(Sph(vec(),python::extract<double>(t[1]))); continue; }
- #if 0
- // compatibility block
- python::extract<python::tuple> tup(t[0]);
- if(tup.check()) { pack.push_back(Sph(Vector3r(python::extract<double>(tup[0]),python::extract<double>(tup[1]),python::extract<double>(tup[2]))),python::extract<double>(t[1])); continue; }
- #endif
- PyErr_SetString(PyExc_TypeError, "List elements must be (Vector3, float)!");
- python::throw_error_already_set();
- }
-};
-
-python::list SpherePack::toList() const {
- python::list ret;
- FOREACH(const Sph& s, pack) ret.append(python::make_tuple(s.c,s.r));
- return ret;
-};
-#if 0
-python::list SpherePack::toList_pointsAsTuples() const {
- python::list ret;
- FOREACH(const Sph& s, pack) ret.append(python::make_tuple(s.c[0],s.c[1],s.c[2],s.r));
- return ret;
-};
-#endif
-
+ const py::tuple& t=py::extract<py::tuple>(l[i]);
+ py::extract<Vector3r> vec(t[0]);
+ if(vec.check()) { pack.push_back(Sph(vec(),py::extract<double>(t[1]),(py::len(t)>2?py::extract<int>(t[2]):-1))); continue; }
+ PyErr_SetString(PyExc_TypeError, "List elements must be (Vector3, float) or (Vector3, float, int)!");
+ py::throw_error_already_set();
+ }
+};
+
+void SpherePack::fromLists(const vector<Vector3r>& centers, const vector<Real>& radii){
+ pack.clear();
+ if(centers.size()!=radii.size()) throw std::invalid_argument(("The same number of centers and radii must be given (is "+lexical_cast<string>(centers.size())+", "+lexical_cast<string>(radii.size())+")").c_str());
+ size_t l=centers.size();
+ for(size_t i=0; i<l; i++){
+ add(centers[i],radii[i]);
+ }
+ cellSize=Vector3r::Zero();
+}
+
+py::list SpherePack::toList() const {
+ py::list ret;
+ FOREACH(const Sph& s, pack) ret.append(s.asTuple());
+ return ret;
+};
void SpherePack::fromFile(string file) {
typedef pair<Vector3r,Real> pairVector3rReal;
@@ -64,7 +62,10 @@
void SpherePack::toFile(const string fname) const {
ofstream f(fname.c_str());
if(!f.good()) throw runtime_error("Unable to open file `"+fname+"'");
- FOREACH(const Sph& s, pack) f<<s.c[0]<<" "<<s.c[1]<<" "<<s.c[2]<<" "<<s.r<<endl;
+ FOREACH(const Sph& s, pack){
+ if(s.clumpId>=0) throw std::invalid_argument("SpherePack with clumps cannot be (currently) exported to a text file.");
+ f<<s.c[0]<<" "<<s.c[1]<<" "<<s.c[2]<<" "<<s.r<<endl;
+ }
f.close();
};
@@ -72,9 +73,10 @@
pack.clear();
Scene* scene=Omega::instance().getScene().get();
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
- shared_ptr<Sphere> intSph=dynamic_pointer_cast<Sphere>(b->shape);
+ if(!b) continue;
+ shared_ptr<Sphere> intSph=dynamic_pointer_cast<Sphere>(b->shape);
if(!intSph) continue;
- pack.push_back(Sph(b->state->pos,intSph->radius));
+ pack.push_back(Sph(b->state->pos,intSph->radius,(b->isClumpMember()?b->clumpId:-1)));
}
if(scene->isPeriodic) { cellSize=scene->cell->getSize(); }
}
@@ -194,8 +196,8 @@
return i;
}
-python::tuple SpherePack::psd(int bins, bool mass) const {
- if(pack.size()==0) return python::make_tuple(python::list(),python::list()); // empty packing
+py::tuple SpherePack::psd(int bins, bool mass) const {
+ if(pack.size()==0) return py::make_tuple(py::list(),py::list()); // empty packing
// find extrema
Real minD=std::numeric_limits<Real>::infinity(); Real maxD=-minD;
// volume, but divided by π*4/3
@@ -211,20 +213,19 @@
if (mass) hist[bin]+=pow(s.r,3)/vol; else hist[bin]+=1./N;
}
for(int i=0; i<bins; i++) cumm[i+1]=min(1.,cumm[i]+hist[i]); // cumm[i+1] is OK, cumm.size()==bins+1
- return python::make_tuple(edges,cumm);
+ return py::make_tuple(edges,cumm);
}
// New code to include the psd giving few points of it
-const float pi = 3.1415926;
long SpherePack::particleSD(Vector3r mn, Vector3r mx, Real rMean, bool periodic, string name, int numSph, const vector<Real>& radii, const vector<Real>& passing, bool passingIsNotPercentageButCount){
vector<Real> numbers;
if(!passingIsNotPercentageButCount){
- Real Vtot=numSph*4./3.*pi*pow(rMean,3.); // total volume of the packing (computed with rMean)
+ Real Vtot=numSph*4./3.*Mathr::PI*pow(rMean,3.); // total volume of the packing (computed with rMean)
// calculate number of spheres necessary per each radius to match the wanted psd
// passing has to contain increasing values
for (size_t i=0; i<radii.size(); i++){
- Real volS=4./3.*pi*pow(radii[i],3.);
+ Real volS=4./3.*Mathr::PI*pow(radii[i],3.);
if (i==0) {numbers.push_back(passing[i]/100.*Vtot/volS);}
else {numbers.push_back((passing[i]-passing[i-1])/100.*Vtot/volS);} //
@@ -273,3 +274,95 @@
return pack.size();
}
+long SpherePack::makeClumpCloud(const Vector3r& mn, const Vector3r& mx, const vector<shared_ptr<SpherePack> >& _clumps, bool periodic, int num){
+ // recenter given clumps and compute their margins
+ vector<SpherePack> clumps; /* vector<Vector3r> margins; */ Vector3r boxMargins(Vector3r::Zero()); Real maxR=0;
+ FOREACH(const shared_ptr<SpherePack>& c, _clumps){
+ SpherePack c2(*c);
+ c2.translate(c2.midPt()); //recenter
+ clumps.push_back(c2);
+ Vector3r cMn,cMx; c2.aabb(cMn,cMx); // centered at zero now, this gives margin
+ //margins.push_back(periodic?cMx:Vector3r::Zero());
+ //boxMargins=boxMargins.cwise().max(cMx);
+ FOREACH(const Sph& s, c2.pack) maxR=max(maxR,s.r); // keep track of maximum sphere radius
+ }
+ Vector3r size=mx-mn;
+ if(periodic)(cellSize=size);
+ const int maxTry=1000;
+ int nGen=0; // number of clumps generated
+ // random point coordinate generator, with non-zero margins if aperiodic
+ static boost::minstd_rand randGen(TimingInfo::getNow(true));
+ typedef boost::variate_generator<boost::minstd_rand&, boost::uniform_real<> > UniRandGen;
+ static UniRandGen rndX(randGen,boost::uniform_real<>(mn[0],mx[0]));
+ static UniRandGen rndY(randGen,boost::uniform_real<>(mn[1],mx[1]));
+ static UniRandGen rndZ(randGen,boost::uniform_real<>(mn[2],mx[2]));
+ static UniRandGen rndUnit(randGen,boost::uniform_real<>(0,1));
+ while(nGen<num || num<0){
+ int clumpChoice=rand()%clumps.size();
+ int tries=0;
+ while(true){ // check for tries at the end
+ Vector3r pos(rndX(),rndY(),rndZ()); // random point
+ // TODO: check this random orientation is homogeneously distributed
+ Quaternionr ori(rndUnit(),rndUnit(),rndUnit(),rndUnit()); ori.normalize();
+ // copy the packing and rotate
+ SpherePack C(clumps[clumpChoice]); C.rotateAroundOrigin(ori); C.translate(pos);
+
+ // find collisions
+ // crude algorithm: check all spheres against all other spheres (slow!!)
+ // use vtkPointLocator, add all existing points and check distance of r+maxRadius, then refine
+ // for periodicity, duplicate points close than boxMargins to the respective boundary
+ if(!periodic){
+ for(size_t i=0; i<C.pack.size(); i++){
+ for(size_t j=0; j<pack.size(); j++){
+ const Vector3r& c(C.pack[i].c); const Real& r(C.pack[i].r);
+ if(pow(r+pack[j].r,2)>=(c-pack[j].c).squaredNorm()) goto overlap;
+ // check that we are not over the box boundary
+ // this could be handled by adjusting the position random interval (by taking off the smallest radius in the clump)
+ // but usually the margin band is relatively small and this does not make the code as hairy
+ if((c+r*Vector3r::Ones()).cwise().max(mx)!=mx || (c-r*Vector3r::Ones()).cwise().min(mn)!=mn) goto overlap;
+ }
+ }
+ }else{
+ for(size_t i=0; i<C.pack.size(); i++){
+ for(size_t j=0; j<pack.size(); j++){
+ const Vector3r& c(C.pack[i].c); const Real& r(C.pack[i].r);
+ Vector3r dr;
+ for(int axis=0; axis<3; axis++) dr[axis]=min(cellWrapRel(c[axis],pack[j].c[axis],pack[j].c[axis]+size[axis]),cellWrapRel(pack[j].c[axis],c[axis],c[axis]+size[axis]));
+ if(pow(pack[j].r+r,2)>= dr.squaredNorm()) goto overlap;
+ }
+ }
+ }
+
+ // add the clump, if no collisions
+ FOREACH(const Sph& s, C.pack){ pack.push_back(Sph(s.c,s.r,/*number clumps consecutively*/nGen)); }
+ nGen++;
+ //cerr<<"O";
+ break; // break away from the try-loop
+
+ overlap:
+ //cerr<<".";
+ if(tries++==maxTry){ // last loop
+ if(num>0) LOG_WARN("Exceeded "<<maxTry<<" attempts to place non-overlapping clump. Only "<<nGen<<" clumps were added, although you requested "<<num);
+ return nGen;
+ }
+ }
+ }
+ return nGen;
+}
+
+bool SpherePack::hasClumps() const { FOREACH(const Sph& s, pack){ if(s.clumpId>=0) return true; } return false; }
+python::tuple SpherePack::getClumps() const{
+ std::map<int,py::list> clumps;
+ py::list standalone; size_t packSize=pack.size();
+ for(size_t i=0; i<packSize; i++){
+ const Sph& s(pack[i]);
+ if(s.clumpId<0) { standalone.append(i); continue; }
+ if(clumps.count(s.clumpId)==0) clumps[s.clumpId]=py::list();
+ clumps[s.clumpId].append(i);
+ }
+ py::list clumpList;
+ typedef std::pair<int,py::list> intListPair;
+ FOREACH(const intListPair& c, clumps) clumpList.append(c.second);
+ return python::make_tuple(standalone,clumpList);
+}
+
=== modified file 'pkg/dem/SpherePack.hpp'
--- pkg/dem/SpherePack.hpp 2010-10-30 17:08:51 +0000
+++ pkg/dem/SpherePack.hpp 2010-11-09 08:27:14 +0000
@@ -19,9 +19,7 @@
#include<yade/lib-base/Logging.hpp>
#include<yade/lib-base/Math.hpp>
-/*! Class representing geometry of spherical packing, with some utility functions.
-
-*/
+/*! Class representing geometry of spherical packing, with some utility functions. */
class SpherePack{
// return coordinate wrapped to x0…x1, relative to x0; don't care about period
// copied from PeriodicInsertionSortCollider
@@ -32,9 +30,13 @@
public:
enum {RDIST_RMEAN, RDIST_POROSITY, RDIST_PSD};
struct Sph{
- Vector3r c; Real r;
- Sph(const Vector3r& _c, Real _r): c(_c), r(_r){};
- python::tuple asTuple() const { return python::make_tuple(c,r); }
+ Vector3r c; Real r; int clumpId;
+ Sph(const Vector3r& _c, Real _r, int _clumpId=-1): c(_c), r(_r), clumpId(_clumpId) {};
+ python::tuple asTuple() const {
+ if(clumpId<0) return python::make_tuple(c,r);
+ return python::make_tuple(c,r,clumpId);
+ }
+ python::tuple asTupleNoClump() const { return python::make_tuple(c,r); }
};
std::vector<Sph> pack;
Vector3r cellSize;
@@ -46,10 +48,8 @@
// I/O
void fromList(const python::list& l);
+ void fromLists(const vector<Vector3r>& centers, const vector<Real>& radii); // used as ctor in python
python::list toList() const;
- #if 0
- python::list toList_pointsAsTuples() const;
- #endif
void fromFile(const string file);
void toFile(const string file) const;
void fromSimulation();
@@ -66,6 +66,10 @@
// interpolate a variable with power distribution (exponent -3) between two margin values, given uniformly distributed x∈(0,1)
Real pow3Interp(Real x,Real a,Real b){ return pow(x*(pow(b,-2)-pow(a,-2))+pow(a,-2),-1./2); }
+ // generate packing of clumps, selected with equal probability
+ // periodic boundary is supported
+ long makeClumpCloud(const Vector3r& mn, const Vector3r& mx, const vector<shared_ptr<SpherePack> >& clumps, bool periodic=false, int num=-1);
+
// periodic repetition
void cellRepeat(Vector3i count);
void cellFill(Vector3r volume);
@@ -85,6 +89,8 @@
return sphVol/(dd[0]*dd[1]*dd[2]);
}
python::tuple psd(int bins=10, bool mass=false) const;
+ bool hasClumps() const;
+ python::tuple getClumps() const;
// transformations
void translate(const Vector3r& shift){ FOREACH(Sph& s, pack) s.c+=shift; }
@@ -92,6 +98,10 @@
if(cellSize!=Vector3r::Zero()) { LOG_WARN("Periodicity reset when rotating periodic packing (non-zero cellSize="<<cellSize<<")"); cellSize=Vector3r::Zero(); }
Vector3r mid=midPt(); Quaternionr q(AngleAxisr(angle,axis)); FOREACH(Sph& s, pack) s.c=q*(s.c-mid)+mid;
}
+ void rotateAroundOrigin(const Quaternionr& rot){
+ if(cellSize!=Vector3r::Zero()){ LOG_WARN("Periodicity reset when rotating periodic packing (non-zero cellSize="<<cellSize<<")"); cellSize=Vector3r::Zero(); }
+ FOREACH(Sph& s, pack) s.c=rot*s.c;
+ }
void scale(Real scale){ Vector3r mid=midPt(); cellSize*=scale; FOREACH(Sph& s, pack) {s.c=scale*(s.c-mid)+mid; s.r*=abs(scale); } }
#if 0
void shrinkMaxRelOverlap(Real maxRelOverlap);
@@ -107,7 +117,7 @@
_iterator iter(){ return *this;}
python::tuple next(){
if(pos==sPack.pack.size()){ PyErr_SetNone(PyExc_StopIteration); python::throw_error_already_set(); }
- return sPack.pack[pos++].asTuple();
+ return sPack.pack[pos++].asTupleNoClump();
}
};
SpherePack::_iterator getIterator() const{ return SpherePack::_iterator(*this);};
=== modified file 'py/pack/_packSpheres.cpp'
--- py/pack/_packSpheres.cpp 2010-09-18 14:10:30 +0000
+++ py/pack/_packSpheres.cpp 2010-11-09 08:27:14 +0000
@@ -6,13 +6,11 @@
BOOST_PYTHON_MODULE(_packSpheres){
python::scope().attr("__doc__")="Creation, manipulation, IO for generic sphere packings.";
YADE_SET_DOCSTRING_OPTS;
- python::class_<SpherePack>("SpherePack","Set of spheres represented as centers and radii. This class is returned by :yref:`yade.pack.randomDensePack`, :yref:`yade.pack.randomPeriPack` and others. The object supports iteration over spheres, as in \n\n\t>>> sp=SpherePack()\n\t>>> for center,radius in sp: print center,radius\n\n\t>>> for sphere in sp: print sphere[0],sphere[1] ## same, but without unpacking the tuple automatically\n\n\t>>> for i in range(0,len(sp)): print sp[i][0], sp[i][1] ## same, but accessing spheres by index\n",python::init<python::optional<python::list> >(python::args("list"),"Empty constructor, optionally taking list [ ((cx,cy,cz),r), … ] for initial data." ))
+ python::class_<SpherePack>("SpherePack","Set of spheres represented as centers and radii. This class is returned by :yref:`yade.pack.randomDensePack`, :yref:`yade.pack.randomPeriPack` and others. The object supports iteration over spheres, as in \n\n\t>>> sp=SpherePack()\n\t>>> for center,radius in sp: print center,radius\n\n\t>>> for sphere in sp: print sphere[0],sphere[1] ## same, but without unpacking the tuple automatically\n\n\t>>> for i in range(0,len(sp)): print sp[i][0], sp[i][1] ## same, but accessing spheres by index\n\n\n.. admonition:: Special constructors\n\n\tConstruct from list of ``[(c1,r1),(c2,r2),…]``. To convert two same-length lists of ``centers`` and ``radii``, construct with ``zip(centers,radii)``.\n",python::init<python::optional<python::list> >(python::args("list"),"Empty constructor, optionally taking list [ ((cx,cy,cz),r), … ] for initial data." ))
.def("add",&SpherePack::add,"Add single sphere to packing, given center as 3-tuple and radius")
.def("toList",&SpherePack::toList,"Return packing data as python list.")
- #if 0
- .def("toList_pointsAsTuples",&SpherePack::toList_pointsAsTuples,"Return packing data as python list, but using only pure-python data types (3-tuples instead of Vector3) (for pickling with cPickle)")
- #endif
.def("fromList",&SpherePack::fromList,"Make packing from given list, same format as for constructor. Discards current data.")
+ .def("fromList",&SpherePack::fromLists,(python::arg("centers"),python::arg("radii")),"Make packing from given list, same format as for constructor. Discards current data.")
.def("load",&SpherePack::fromFile,(python::arg("fileName")),"Load packing from external text file (current data will be discarded).")
.def("save",&SpherePack::toFile,(python::arg("fileName")),"Save packing to external text file (will be overwritten).")
.def("fromSimulation",&SpherePack::fromSimulation,"Make packing corresponding to the current simulation. Discards current data.")
@@ -22,6 +20,7 @@
.def("psd",&SpherePack::psd,(python::arg("bins")=10,python::arg("mass")=false),"Return `particle size distribution <http://en.wikipedia.org/wiki/Particle_size_distribution>`__ of the packing.\n:param int bins: number of bins between minimum and maximum diameter\n:param mass: Compute relative mass rather than relative particle count for each bin. Corresponds to :yref:`distributeMass parameter for makeCloud<yade.pack.SpherePack.makeCloud>`.\n:returns: tuple of ``(cumm,edges)``, where ``cumm`` are cummulative fractions for respective diameters and ``edges`` are those diameter values. Dimension of both arrays is equal to ``bins+1``.")
// new psd
.def("particleSD",&SpherePack::particleSD,(python::arg("minCorner"),python::arg("maxCorner"),python::arg("rMean"),python::arg("periodic")=false,python::arg("name"),python::arg("numSph"),python::arg("radii")=vector<Real>(),python::arg("passing")=vector<Real>(),python::arg("passingIsNotPercentageButCount")=false),"Create random packing enclosed in box given by minCorner and maxCorner, containing numSph spheres. Returns number of created spheres, which can be < num if the packing is too tight. The computation is done according to the given psd.")
+ .def("makeClumpCloud",&SpherePack::makeClumpCloud,(python::arg("minCorner"),python::arg("maxCorner"),python::arg("clumps"),python::arg("periodic")=false,python::arg("num")=-1),"Create random loose packing of clumps within box given by *minCorner* and *maxCorner*. Clumps are selected with equal probability. At most *num* clumps will be positioned if *num* is positive; otherwise, as many clumps as possible will be put in space, until maximum number of attemps to place a new clump randomly is attained.")
//
.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.")
@@ -33,6 +32,8 @@
.def("translate",&SpherePack::translate,"Translate all spheres by given vector.")
.def("rotate",&SpherePack::rotate,(python::arg("axis"),python::arg("angle")),"Rotate all spheres around packing center (in terms of aabb()), given axis and angle of the rotation.")
.def("scale",&SpherePack::scale,"Scale the packing around its center (in terms of aabb()) by given factor (may be negative).")
+ .def("hasClumps",&SpherePack::hasClumps,"Whether this object contains clumps.")
+ .def("getClumps",&SpherePack::getClumps,"Return lists of sphere ids sorted by clumps they belong to. The return value is (standalones,[clump1,clump2,…]), where each item is list of id's of spheres.")
.def("__len__",&SpherePack::len,"Get number of spheres in the packing")
.def("__getitem__",&SpherePack::getitem,"Get entry at given index, as tuple of center and radius.")
.def("__iter__",&SpherePack::getIterator,"Return iterator over spheres.")
=== modified file 'py/pack/pack.py'
--- py/pack/pack.py 2010-10-29 10:12:44 +0000
+++ py/pack/pack.py 2010-11-09 08:27:14 +0000
@@ -92,7 +92,19 @@
O.periodic=True
O.cell.trsf=rot
O.cell.refSize=self.cellSize
- return O.bodies.append([utils.sphere(rot*c,r,**kw) for c,r in self])
+ 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(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[i].shape.color=O.bodies[clump[0]].shape.color
+ return ids+clumpIds
SpherePack.toSimulation=SpherePack_toSimulation
=== modified file 'py/wrapper/customConverters.cpp'
--- py/wrapper/customConverters.cpp 2010-10-29 10:12:44 +0000
+++ py/wrapper/customConverters.cpp 2010-11-09 08:27:14 +0000
@@ -35,6 +35,7 @@
#include<yade/pkg-common/Dispatching.hpp>
#include<yade/pkg-common/Callbacks.hpp>
+#include<yade/pkg-dem/SpherePack.hpp>
#ifdef YADE_OPENGL
#include<yade/pkg-common/GLDrawFunctors.hpp>
#include<yade/pkg-common/OpenGLRenderer.hpp>
@@ -225,6 +226,7 @@
VECTOR_SEQ_CONV(shared_ptr<LawFunctor>);
VECTOR_SEQ_CONV(shared_ptr<IntrCallback>);
VECTOR_SEQ_CONV(shared_ptr<BodyCallback>);
+ VECTOR_SEQ_CONV(shared_ptr<SpherePack>);
#ifdef YADE_OPENGL
VECTOR_SEQ_CONV(shared_ptr<GlBoundFunctor>);
VECTOR_SEQ_CONV(shared_ptr<GlStateFunctor>);
=== added file 'scripts/test/clumpPack.py'
--- scripts/test/clumpPack.py 1970-01-01 00:00:00 +0000
+++ scripts/test/clumpPack.py 2010-11-09 08:27:14 +0000
@@ -0,0 +1,7 @@
+# create a few clump configurations by hand
+from yade import pack
+c1=pack.SpherePack([((0,0,0),.5),((.5,0,0),.5),((0,.5,0),.3)])
+c2=pack.SpherePack([((0,0,0),.5),((.7,0,0),.3),((.9,0,0),.2)])
+sp=pack.SpherePack()
+sp.makeClumpCloud((0,0,0),(10,10,10),[c1,c2],periodic=False)
+sp.toSimulation()
