yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #03904
[Branch ~yade-dev/yade/trunk] Rev 2136: 1. Add toy grid-based collider (very slow)
------------------------------------------------------------
revno: 2136
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Fri 2010-04-09 13:26:20 +0200
message:
1. Add toy grid-based collider (very slow)
added:
pkg/dem/Engine/GlobalEngine/FlatGridCollider.cpp
pkg/dem/Engine/GlobalEngine/FlatGridCollider.hpp
scripts/test/flat-collider.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
=== added file 'pkg/dem/Engine/GlobalEngine/FlatGridCollider.cpp'
--- pkg/dem/Engine/GlobalEngine/FlatGridCollider.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Engine/GlobalEngine/FlatGridCollider.cpp 2010-04-09 11:26:20 +0000
@@ -0,0 +1,105 @@
+// 2010 © Václav Šmilauer <eudoxos@xxxxxxxx>
+#include<yade/pkg-dem/FlatGridCollider.hpp>
+#include<yade/core/Scene.hpp>
+
+#include<yade/pkg-common/Sphere.hpp>
+#include<yade/pkg-dem/NewtonIntegrator.hpp>
+//#include<yade/pkg-common/Facet.hpp>
+
+YADE_PLUGIN((FlatGridCollider));
+CREATE_LOGGER(FlatGridCollider);
+
+bool FlatGridCollider::isActivated(){
+ // keep interactions trequested for deletion as potential (forget removal requests)
+ scene->interactions->clearPendingErase();
+ if(!newton) return true;
+ // handle verlet distance
+ fastestBodyMaxDist+=sqrt(newton->maxVelocitySq)*scene->dt;
+ if(fastestBodyMaxDist>verletDist) return true;
+ return false;
+}
+
+void FlatGridCollider::action(){
+ if(!newton){
+ FOREACH(const shared_ptr<Engine>& e, scene->engines){ newton=dynamic_pointer_cast<NewtonIntegrator>(e); if(newton) break; }
+ if(!newton){ throw runtime_error("FlatGridCollider: Unable to find NewtonIntegrator in engines."); }
+ }
+ fastestBodyMaxDist=0;
+ // make interaction loop delete unseen potential interactions
+ scene->interactions->iterColliderLastRun=scene->currentIteration;
+ // adjust grid if necessary
+ updateGrid();
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if(!b) continue; // deleted bodies
+ updateBodyCells(b);
+ }
+ updateCollisions();
+}
+
+
+// called from the ctor
+void FlatGridCollider::initIndices(){
+ sphereIdx=facetIdx=wallIdx=boxIdx=-1;
+ sphereIdx=Sphere::getClassIndexStatic();
+ LOG_DEBUG("sphereIdx="<<sphereIdx);
+}
+
+void FlatGridCollider::updateGrid(){
+ // checks
+ if(step<=0) throw std::runtime_error("FlatGridCollider::step must be positive.");
+ Vector3r aabbSize=aabbMax-aabbMin; if(aabbSize[0]<=0 || aabbSize[1]<=0 || aabbSize[2]<=0) throw std::runtime_error("FlatGridCollider::{aabbMin,aabbMax} must give positive volume.");
+ // compute new size
+ grid.step=step;
+ grid.mn=aabbMin;
+ for(int i=0;i<3;i++)grid.size[i]=(int)ceil((aabbMax[i]-aabbMin[i])/step);
+ grid.mx=aabbMin+Vector3r(grid.size[0]*step,grid.size[1]*step,grid.size[2]*step);
+ size_t len=grid.size[0]*grid.size[1]*grid.size[2];
+ // reset grid data
+ grid.data.clear(); grid.data.resize(len); // delete and recreate; could be optimized somehow
+ LOG_DEBUG("New grid "<<grid.size[0]<<"×"<<grid.size[1]<<"×"<<grid.size[2]<<"="<<len<<" cells, step "<<step<<"; spans ("<<grid.mn<<")--("<<grid.mx<<").");
+}
+
+void FlatGridCollider::updateBodyCells(const shared_ptr<Body>& b){
+ if(!b->shape) return; const Shape* shape(b->shape.get());
+ // Sphere
+ if(shape->getClassIndex()==sphereIdx){
+ Real r=static_cast<const Sphere*>(shape)->radius+verletDist;
+ const Vector3r& C=b->state->pos;
+ // create "bounding cells" and traverse them one by one; integrized coords can be _outside_ grid, they will be forced to closest cells before insertion
+ Vector3i cMn=grid.pt2int(C-Vector3r(r,r,r)), cMx=grid.pt2int(C+Vector3r(r,r,r)), cC=grid.pt2int(C), cPt;
+ LOG_TRACE("Sphere #"<<b->id<<", C="<<C<<", cMn="<<cMn<<", cC="<<cC<<", cMx="<<cMx);
+ for(cPt[0]=cMn[0]; cPt[0]<=cMx[0]; cPt[0]++) for(cPt[1]=cMn[1]; cPt[1]<=cMx[1]; cPt[1]++) for(cPt[2]=cMn[2]; cPt[2]<=cMx[2]; cPt[2]++){
+ // find closest cell point (to cC); keep coordinate in same line (cPt[i]==cC[i]); take upper/lower point in lower/upper lines (cPt[i]<cC[i])
+ Vector3r ccp;
+ for(int i=0;i<3;i++) ccp[i]=(cPt[i]==cC[i] ? C[i] : (grid.mn[i]+grid.step*(cPt[i]+(cPt[i]<cC[i] ? 1 : 0))));
+ if((C-ccp).SquaredLength()<=r*r){ // closest cell point it inside the spehre; add the sphere to cell at cell position
+ Vector3i cPtIn(grid.fitGrid(cPt));
+ // perhaps slower, but inserts each body only once into the cell (meaningful if outside grid: multiple integer coords coolapse in one cell)
+ { Grid::idVector& vv=grid(cPtIn); if(vv.size()==0 || *(vv.rbegin())!=b->id) vv.push_back(b->id);}
+ //grid(cPtIn).push_back(b->id);
+ LOG_TRACE("Added sphere #"<<b->id<<" to cell ("<<cPtIn<<")â†["<<cPt<<"]; center ("<<C<<"), closest ("<<ccp<<"), dist "<<(C-ccp).Length());
+ }
+ }
+ return;
+ }
+ throw std::runtime_error("FlatGridCollider::updateBodyCells does not handle Shape of type "+shape->getClassName()+"!");
+}
+
+void FlatGridCollider::updateCollisions(){
+ shared_ptr<InteractionContainer>& intrs=scene->interactions;
+ const long& iter=scene->currentIteration;
+ // create interactions for all combinations of bodies within one cell
+ FOREACH(const Grid::idVector& v, grid.data){
+ size_t sz=v.size();
+ for(size_t i=0; i<sz; i++) for(size_t j=i+1; j<sz; j++){
+ body_id_t id1=v[i], id2=v[j];
+ if(id1==id2) continue; // at grid boundary, it is possible to have one body more times in one cell
+ const shared_ptr<Interaction>& I=intrs->find(id1,id2);
+ if(I){ I->iterLastSeen=iter; continue; }
+ // no interaction yet
+ if(!Collider::mayCollide(Body::byId(id1,scene).get(),Body::byId(id2,scene).get())) continue;
+ intrs->insert(shared_ptr<Interaction>(new Interaction(id1,id2)));
+ LOG_TRACE("Created new interaction #"<<id1<<"+#"<<id2);
+ }
+ }
+}
=== added file 'pkg/dem/Engine/GlobalEngine/FlatGridCollider.hpp'
--- pkg/dem/Engine/GlobalEngine/FlatGridCollider.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Engine/GlobalEngine/FlatGridCollider.hpp 2010-04-09 11:26:20 +0000
@@ -0,0 +1,46 @@
+// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
+#include<vector>
+#include<yade/core/Collider.hpp>
+class NewtonIntegrator;
+class FlatGridCollider: public Collider{
+ struct Grid{
+ typedef std::vector<body_id_t> idVector;
+ Vector3i size;
+ Vector3r mn, mx;
+ Real step;
+ // convert point into its integral coordinates (can be outside grid, use fitGrid to coords inside)
+ Vector3i pt2int(const Vector3r& pt){ Vector3i ret; for(int i=0;i<3;i++)ret[i]=floor((pt[i]-mn[1])/step); return ret; }
+ std::vector<idVector> data;
+ // force integral coordinate inside (0,sz-1)
+ int fit(int i, int sz) const { return max(0,min(i,sz-1)); }
+ Vector3i fitGrid(const Vector3i& v){ return Vector3i(fit(v[0],size[0]),fit(v[1],size[1]),fit(v[2],size[2])); }
+ // linearize x,y,z → n in data vector
+ size_t lin(int x, int y, int z) const { return fit(x,size[0])+size[0]*fit(y,size[1])+size[0]*size[1]*fit(z,size[2]); }
+ // return vector of ids at (x,y,z)
+ idVector& operator()(int x, int y, int z){ return data[lin(x,y,z)];}
+ idVector& operator()(const Vector3i& v){ return data[lin(v[0],v[1],v[2])];}
+ const idVector& operator()(int x, int y, int z) const { return data[lin(x,y,z)];}
+ };
+ Grid grid;
+ int sphereIdx, facetIdx, wallIdx, boxIdx;
+ // needed for maxVelSq
+ shared_ptr<NewtonIntegrator> newton;
+ // track maximum distance of fastest body, at every step in isActivated
+ Real fastestBodyMaxDist;
+ void initIndices();
+ void updateGrid();
+ void updateBodyCells(const shared_ptr<Body>& b);
+ void updateCollisions();
+ virtual void action();
+ virtual bool isActivated();
+ DECLARE_LOGGER;
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(FlatGridCollider,Collider,"Non-optimized grid collider, storing grid as dense flat array. Each body is assigned to (possibly multiple) cells, which are arranged in regular grid between *aabbMin* and *aabbMax*, with cell size *step* (same in all directions). Bodies outsize (*aabbMin*, *aabbMax*) are handled gracefully, assigned to closest cells (this will create spurious potential interactions). *verletDist* determines how much is each body enlarged to avoid collision detection at every step.\n\n.. note::\n\tThis collider keeps all cells in linear memory array, therefore will be memory-inefficient for sparse simulations.\n\n.. warning::\n\t:yref:`Body::bound` objects are not used, :yref:`BoundFunctors<BoundFunctor>` are not used either: assigning cells to bodies is hard-coded internally. Currently handles :yref:`Shapes<Shape>` are: :yref:`Sphere`.\n\n.. note::\n\tPeriodic boundary is not handled (yet).\n\n",
+ ((Real,verletDist,0,"Length by which enlarge space occupied by each particle; avoids running collision detection at every step."))
+ ((Vector3r,aabbMin,Vector3r::ZERO,"Lower corner of grid."))
+ ((Vector3r,aabbMax,Vector3r::ZERO,"Upper corner of grid (approximate, might be rouded up to *minStep*."))
+ ((Real,step,0,"Step in the grid (cell size)")),
+ initIndices();
+ );
+};
+REGISTER_SERIALIZABLE(FlatGridCollider);
=== added file 'scripts/test/flat-collider.py'
--- scripts/test/flat-collider.py 1970-01-01 00:00:00 +0000
+++ scripts/test/flat-collider.py 2010-04-09 11:26:20 +0000
@@ -0,0 +1,28 @@
+from yade import log,utils,pack,timing
+#log.setLevel('FlatGridCollider',log.TRACE)
+#O.bodies.append([ utils.sphere((0.2,0,0),.5,dynamic=False), utils.sphere((0.2,0.0,1.01),.5), ])
+O.bodies.append(pack.regularHexa(pack.inAlignedBox((0,0,0),(10,10,1)),radius=.5,gap=0,dynamic=False))
+O.bodies.append(pack.regularOrtho(pack.inAlignedBox((3,3,3),(7,7,4)),radius=.05,gap=0))
+O.engines=[
+ ForceResetter(),
+ FlatGridCollider(step=.2,aabbMin=(0,0,0),aabbMax=(10,10,5),verletDist=0.005),
+ #BoundDispatcher([Bo1_Sphere_Aabb()]), InsertionSortCollider(sweepLength=0.005),
+ InteractionDispatchers(
+ [Ig2_Sphere_Sphere_Dem3DofGeom()],
+ [Ip2_FrictMat_FrictMat_FrictPhys()],
+ [Law2_Dem3DofGeom_FrictPhys_Basic()],
+ ),
+ GravityEngine(gravity=[0,0,-10]),
+ NewtonIntegrator(damping=0.4),
+]
+O.dt=.6*utils.PWaveTimeStep()
+O.saveTmp()
+#O.step()
+#while True:
+# O.step()
+# if len(O.interactions)>0 or O.bodies[1].state.pos[2]<.97: break
+O.timingEnabled=True
+O.run(5000,True)
+timing.stats()
+import sys
+#sys.exit(0)
Follow ups
