yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06125
[Branch ~yade-dev/yade/trunk] Rev 2536: merge with trunk
Merge authors:
jduriez <jduriez@c1solimara-l>
Václav Šmilauer (eudoxos)
------------------------------------------------------------
revno: 2536 [merge]
committer: Sergei D. <sega@think>
branch nick: trunk
timestamp: Tue 2010-11-09 13:58:43 +0300
message:
merge with trunk
added:
scripts/test/clumpPack.py
modified:
doc/yade-articles.bib
pkg/common/OpenGLRenderer.cpp
pkg/common/OpenGLRenderer.hpp
pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp
pkg/dem/NormalInelasticityLaw.cpp
pkg/dem/NormalInelasticityLaw.hpp
pkg/dem/NormalInelasticityPhys.hpp
pkg/dem/SpherePack.cpp
pkg/dem/SpherePack.hpp
py/pack/_packSpheres.cpp
py/pack/pack.py
py/wrapper/customConverters.cpp
scripts/normalInelasticity-test.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 '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/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp'
--- pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp 2010-11-08 15:26:46 +0000
@@ -67,18 +67,9 @@
contactPhysics->prevNormal = geom->normal;
contactPhysics->knLower = Kn;
- contactPhysics->kn = Kn;
-
+ contactPhysics->kn = Kn;
contactPhysics->ks = Ks;
- contactPhysics->initialOrientation1 = Body::byId(interaction->getId1())->state->ori;
- contactPhysics->initialOrientation2 = Body::byId(interaction->getId2())->state->ori;
- contactPhysics->initialPosition1 = Body::byId(interaction->getId1())->state->pos;
- contactPhysics->initialPosition2 = Body::byId(interaction->getId2())->state->pos;
contactPhysics->kr = Kr;
- contactPhysics->initialContactOrientation.setFromTwoVectors(Vector3r(1.0,0.0,0.0),geom->normal);
- contactPhysics->currentContactOrientation = contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact1 = contactPhysics->initialOrientation1.conjugate() * contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact2 = contactPhysics->initialOrientation2.conjugate() * contactPhysics->initialContactOrientation;
}
}
=== modified file 'pkg/dem/NormalInelasticityLaw.cpp'
--- pkg/dem/NormalInelasticityLaw.cpp 2010-11-03 15:50:02 +0000
+++ pkg/dem/NormalInelasticityLaw.cpp 2010-11-08 15:26:46 +0000
@@ -165,28 +165,10 @@
if(momentRotationLaw)
{
-// {// updates only orientation of contact (local coordinate system)
-// Vector3r axis = currentContactPhysics->prevNormal.cross(currentContactGeometry->normal); axis.normalize();
-// Real angle = unitVectorsAngle(currentContactPhysics->prevNormal,currentContactGeometry->normal);
-// Quaternionr align(AngleAxisr(angle,axis));
-// currentContactPhysics->currentContactOrientation = align * currentContactPhysics->currentContactOrientation;
-// }
-//
-// Quaternionr delta( de1->se3.orientation * currentContactPhysics->initialOrientation1.conjugate() *
-// currentContactPhysics->initialOrientation2 * de2->se3.orientation.conjugate());
-//
-// AngleAxisr aa(delta); // aa.axis() of rotation - this is the Moment direction UNIT vector; angle represents the power of resistant ELASTIC moment
-// if(angle > Mathr::PI) angle -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
-
-//Real elasticMoment = currentContactPhysics->kr * std::abs(angle); // positive value (*)
-
-// Real angle_twist(aa.angle() * aa.axis().dot(currentContactGeometry->normal) );
-// Vector3r axis_twist(angle_twist * currentContactGeometry->normal);
- Vector3r moment_twist( (currentContactGeometry->getTwist()*currentContactPhysics->kr)*currentContactGeometry->normal );
-// Vector3r axis_bending(aa.angle()*aa.axis() - axis_twist);
- Vector3r moment_bending( currentContactGeometry->getBending() * currentContactPhysics->kr );
-
- Vector3r moment = moment_twist + moment_bending;
+ currentContactPhysics->moment_twist = (currentContactGeometry->getTwist()*currentContactPhysics->kr)*currentContactGeometry->normal ;
+ currentContactPhysics->moment_bending = currentContactGeometry->getBending() * currentContactPhysics->kr;
+
+ moment = currentContactPhysics->moment_twist + currentContactPhysics->moment_bending;
// Limitation by plastic threshold
if (!momentAlwaysElastic)
@@ -196,7 +178,6 @@
if(normeMoment>normeMomentMax)
{
moment *= normeMomentMax/normeMoment;
-// nbreInteracMomPlastif++;
}
}
scene->forces.addTorque(id1,-moment);
@@ -205,8 +186,6 @@
// ******** Moment law END ******* //
currentContactPhysics->prevNormal = currentContactGeometry->normal;
-// }
-// }
}
}
=== modified file 'pkg/dem/NormalInelasticityLaw.hpp'
--- pkg/dem/NormalInelasticityLaw.hpp 2010-11-03 15:50:02 +0000
+++ pkg/dem/NormalInelasticityLaw.hpp 2010-11-08 15:26:46 +0000
@@ -18,16 +18,20 @@
class Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity : public LawFunctor
{
+ private :
+ Vector3r moment; // contact torque of the interaction
public :
virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
FUNCTOR2D(ScGeom,NormalInelasticityPhys);
- YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity,
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity,
LawFunctor,
- "Contact law used to simulate granulate filler in rock joints [Duriez2009a]_, [Duriez2010]_. It includes possibility of cohesion, moment transfer and inelastic compression behaviour (to reproduce the normal inelasticity observed for rock joints, for the latter).\n\n The moment transfer relation corresponds to the adaptation of the work of Plassiard & Belheine (see in [DeghmReport2006]_ for example), which was realized by J. Kozicki, and is now coded in ScGeom6D.\n\n As others :yref:`LawFunctor`, it uses pre-computed data of the interactions (rigidities, friction angles -with their tan()-, orientations of the interactions); this work is done here in :yref:`Ip2_2xNormalInelasticMat_NormalInelasticityPhys`.\n\n To use this you should also use :yref:`NormalInelasticMat` as material type of the bodies.\n\n The effects of this law are illustrated in scripts/normalInelasticityTest.py",
+ "Contact law used to simulate granulate filler in rock joints [Duriez2009a]_, [Duriez2010]_. It includes possibility of cohesion, moment transfer and inelastic compression behaviour (to reproduce the normal inelasticity observed for rock joints, for the latter).\n\n The moment transfer relation corresponds to the adaptation of the work of Plassiard & Belheine (see in [DeghmReport2006]_ for example), which was realized by J. Kozicki, and is now coded in :yref:`ScGeom6D`.\n\n As others :yref:`LawFunctor`, it uses pre-computed data of the interactions (rigidities, friction angles -with their tan()-, orientations of the interactions); this work is done here in :yref:`Ip2_2xNormalInelasticMat_NormalInelasticityPhys`.\n\n To use this you should also use :yref:`NormalInelasticMat` as material type of the bodies.\n\n The effects of this law are illustrated in scripts/normalInelasticityTest.py",
((bool,momentRotationLaw,true,,"boolean, true=> computation of a torque (against relative rotation) exchanged between particles"))
((bool,momentAlwaysElastic,false,,"boolean, true=> the torque (computed only if momentRotationLaw !!) is not limited by a plastic threshold"))
+ ,
+ moment=Vector3r::Zero();
);
};
=== modified file 'pkg/dem/NormalInelasticityPhys.hpp'
--- pkg/dem/NormalInelasticityPhys.hpp 2010-11-03 15:50:02 +0000
+++ pkg/dem/NormalInelasticityPhys.hpp 2010-11-08 15:26:46 +0000
@@ -17,23 +17,19 @@
virtual ~NormalInelasticityPhys();
YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticityPhys,FrictPhys,
- "Physics (of interaction) for using :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity` : with inelastic unloadings",
- ((Real,unMax,0.0,,"the maximum value of penetration depth of the history of this interaction"))
- ((Real,previousun,0.0,,"the value of this un at the last time step"))
- ((Real,previousFn,0.0,,"the value of the normal force at the last time step"))
- ((Quaternionr,initialOrientation1,Quaternionr::Identity(),,""))
- ((Quaternionr,initialOrientation2,Quaternionr::Identity(),,""))
- ((Quaternionr,orientationToContact1,Quaternionr::Identity(),,""))
- ((Quaternionr,orientationToContact2,Quaternionr::Identity(),,""))
- ((Quaternionr,currentContactOrientation,Quaternionr::Identity(),,""))
- ((Quaternionr,initialContactOrientation,Quaternionr::Identity(),,""))
- ((Vector3r,initialPosition1,Vector3r::Zero(),,""))
- ((Vector3r,initialPosition2,Vector3r::Zero(),,""))
- ((Real,forMaxMoment,1.0,,"parameter stored for each interaction, and allowing to compute the maximum value of the exchanged torque : TorqueMax= forMaxMoment * NormalForce"))
- ((Real,kr,0.0,,"the rolling stiffness of the interaction"))
- ((Real,knLower,0.0,,"the stifness corresponding to a virgin load for example")),
- createIndex();
- );
+ "Physics (of interaction) for using :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity` : with inelastic unloadings",
+ ((Real,unMax,0.0,,"the maximum value of penetration depth of the history of this interaction"))
+ ((Real,previousun,0.0,,"the value of this un at the last time step"))
+ ((Real,previousFn,0.0,,"the value of the normal force at the last time step"))
+ ((Real,forMaxMoment,1.0,,"parameter stored for each interaction, and allowing to compute the maximum value of the exchanged torque : TorqueMax= forMaxMoment * NormalForce"))
+ ((Real,kr,0.0,,"the rolling stiffness of the interaction"))
+ ((Real,knLower,0.0,,"the stifness corresponding to a virgin load for example"))
+ // internal attributes
+ ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Twist moment. Defined here, being initialized as it should be, to be used in :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`"))
+ ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Bending moment. Defined here, being initialized as it should be, to be used in :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`"))
+ ,
+ createIndex();
+ );
REGISTER_CLASS_INDEX(NormalInelasticityPhys,FrictPhys);
};
=== 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>);
=== modified file 'scripts/normalInelasticity-test.py'
--- scripts/normalInelasticity-test.py 2010-11-03 15:50:02 +0000
+++ scripts/normalInelasticity-test.py 2010-11-08 15:26:46 +0000
@@ -28,7 +28,7 @@
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
- [Ig2_Sphere_Sphere_ScGeom()],
+ [Ig2_Sphere_Sphere_ScGeom6D()],
[Ip2_2xNormalInelasticMat_NormalInelasticityPhys(betaR=0.24)],
[Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity()]
),
@@ -52,7 +52,7 @@
vecFn=i.phys.normalForce
vecDist=upperSphere.state.pos-lowerSphere.state.pos
plot.addData(normFn=vecFn.norm(),normFnBis=vecFn.norm(),fnY=vecFn[1],step=O.iter,
- unPerso=lowerSphere.shape.radius+upperSphere.shape.radius-vecDist.norm(),unVrai=i.geom.penetrationDepth,
+ unPerso=lowerSphere.shape.radius+upperSphere.shape.radius-vecDist.norm(),unTrue=i.geom.penetrationDepth,
gamma=upperSphere.state.pos[0]-lowerSphere.state.pos[0],fx=O.forces.f(0)[0],torque=O.forces.t(1)[2])
#print i.geom.penetrationDepth
@@ -65,17 +65,17 @@
O.run(2,True) #cycles "for free", so that the interaction between spheres will be defined (with his physics and so on)
O.engines=O.engines+[PyRunner(iterPeriod=1,command='defData()')]
-
-
O.run(40,True)
+print 'End of normal loading'
+
# define of the plots to be made : un(step), and Fn(un)
-plot.plots={'step':('unVrai',),'unPerso':('normFn',),'unVrai':('normFnBis',)}
+plot.plots={'step':('unTrue','torque',),'unPerso':('normFn',),'unTrue':('normFnBis',)}
plot.plot()
raw_input()
-print 'End of normal loading'
+print 'Type Return to go ahead'
print ''
-#NB : these different unVrai and unPerso illustrate the definition of penetrationDepth really used in the code (computed in Ig2_Sphere_Sphere_ScGeom) which is slightly different from R1 + R2 - Distance (see for example this "shift2"). According to the really used penetrationDepth, Fn evolves as it should
+#NB : these different unTrue and unPerso illustrate the definition of penetrationDepth really used in the code (computed in Ig2_Sphere_Sphere_ScGeom) which is slightly different from R1 + R2 - Distance (see for example this "shift2"). According to the really used penetrationDepth, Fn evolves as it should
#O.saveTmp('EndComp')
#O.save('FinN_I_Test.xml')
@@ -89,15 +89,17 @@
Vector3.__init__(dpos,1*O.dt,0,0)
O.engines=O.engines[:3]+[StepDisplacer(ids=[1],mov=dpos,setVelocities=True)]+O.engines[4:]
O.run(1000)
-plot.plots={'step':('gamma',),'gamma':('fx',)}
+print 'End of tangential loading'
+plot.plots={'step':('gamma','torque',),'gamma':('fx',)}
plot.plot()
raw_input()
+print 'Type Return to go ahead'
plot.plots={'normFn':('fx',)}
plot.plot()
raw_input()
+print 'Type Return to go ahead'
#pylab.show() #to pause on the plot window. Effective only first time
-print 'End of tangential loading'
print ''
#-- Comments (r2528) --#
@@ -117,7 +119,8 @@
upperSphere.state.angVel=Vector3(0,0,1)
upperSphere.state.vel=Vector3(0,0,0)
i=O.interactions[1,0]
-O.engines=O.engines[:4]+[NewtonIntegrator()]+O.engines[5:]#+[PyRunner(iterPeriod=1,command='printInfo()')]
+
+O.engines=O.engines[:3]+[NewtonIntegrator()]+O.engines[4:]#+[PyRunner(iterPeriod=1,command='printInfo()')]
def printInfo():
@@ -126,7 +129,10 @@
print i.geom.penetrationDepth
O.run(8000,True)
+print 'End of rotationnal loading'
+
plot.plots={'step':('torque',)}
plot.plot()
+print 'Type Return to go ahead'
#-- Comments : TO DO
=== 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()
