yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01623
[svn] r1928 - in trunk: core extra pkg/common/Engine/StandAloneEngine py/yadeWrapper scripts/test
Author: eudoxos
Date: 2009-08-07 00:05:28 +0200 (Fri, 07 Aug 2009)
New Revision: 1928
Added:
trunk/extra/PeriodicInsertionSortCollider.cpp
trunk/extra/PeriodicInsertionSortCollider.hpp
trunk/scripts/test/periodic-simple.py
Modified:
trunk/core/Interaction.cpp
trunk/core/Interaction.hpp
trunk/core/InteractionContainer.hpp
trunk/core/MetaBody.cpp
trunk/core/MetaBody.hpp
trunk/core/yade.cpp
trunk/extra/SConscript
trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp
trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp
trunk/py/yadeWrapper/yadeWrapper.cpp
Log:
1. Infrastructure for periodic collision handling. Not at all functional at this moment, don't even try to use.
Modified: trunk/core/Interaction.cpp
===================================================================
--- trunk/core/Interaction.cpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/core/Interaction.cpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -21,6 +21,7 @@
isNeighbor = true;//NOTE : TriangulationCollider needs that
iterMadeReal=-1;
functorCache.geomExists=true;
+ cellDist=Vector3<int>(0,0,0);
//functorCache.geom=shared_ptr<InteractionGeometryEngineUnit>(); functorCache.phys=shared_ptr<InteractionPhysicsEngineUnit>(); functorCache.constLaw=shared_ptr<ConstitutiveLaw>();
}
Modified: trunk/core/Interaction.hpp
===================================================================
--- trunk/core/Interaction.hpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/core/Interaction.hpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -38,6 +38,10 @@
//! NOTE : TriangulationCollider needs this (nothing else)
bool isNeighbor;
+ //! relative distance between bodies, given in (MetaBody::cellMax-MetaBody::cellMin) units
+ //! Position of id1 must be incremented by that distance so that there is spatial interaction
+ Vector3<int> cellDist;
+
shared_ptr<InteractionGeometry> interactionGeometry;
shared_ptr<InteractionPhysics> interactionPhysics;
@@ -74,6 +78,7 @@
(iterMadeReal)
(interactionGeometry)
(interactionPhysics)
+ (cellDist)
);
REGISTER_CLASS_AND_BASE(Interaction,Serializable);
};
Modified: trunk/core/InteractionContainer.hpp
===================================================================
--- trunk/core/InteractionContainer.hpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/core/InteractionContainer.hpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -146,13 +146,13 @@
Returns number of interactions, have they been erased or not (this is useful to check if there were some erased, after traversing those)
*/
- template<class T> int erasePending(const T& t){
+ template<class T> int erasePending(const T& t, MetaBody* rb){
int ret=0;
#ifdef YADE_OPENMP
// shadow the this->pendingErase by the local variable, to share the code
FOREACH(list<bodyIdPair>& pendingErase, threadsPendingErase){
#endif
- FOREACH(const Vector2<body_id_t>& p, pendingErase){ ret++; if(t.shouldBeErased(p[0],p[1])) erase(p[0],p[1]); }
+ FOREACH(const Vector2<body_id_t>& p, pendingErase){ ret++; if(t.shouldBeErased(p[0],p[1],rb)) erase(p[0],p[1]); }
pendingErase.clear();
#ifdef YADE_OPENMP
}
Modified: trunk/core/MetaBody.cpp
===================================================================
--- trunk/core/MetaBody.cpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/core/MetaBody.cpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -47,6 +47,7 @@
isDynamic=false;
dt=1e-8;
selectedBody=-1;
+ isPeriodic=false;
// FIXME: move MetaInteractingGeometry to core and create it here right away
// interactingGeometry=shared_ptr<InteractingGeometry>(new MetaInteractingGeometry);
Modified: trunk/core/MetaBody.hpp
===================================================================
--- trunk/core/MetaBody.hpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/core/MetaBody.hpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -51,6 +51,8 @@
long stopAtIteration;
Real stopAtVirtTime;
Real stopAtRealTime;
+ Vector3r cellMin, cellMax;
+ bool isPeriodic;
bool needsInitializers;
// for GL selection
@@ -70,6 +72,9 @@
(currentIteration)
(simulationTime)
(stopAtIteration)
+ (cellMin)
+ (cellMax)
+ (isPeriodic)
);
REGISTER_CLASS_AND_BASE(MetaBody,Body);
};
Modified: trunk/core/yade.cpp
===================================================================
--- trunk/core/yade.cpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/core/yade.cpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -262,6 +262,7 @@
LOG_INFO("Loading plugins");
Omega::instance().scanPlugins(string(PREFIX "/lib/yade" SUFFIX "/plugins" ));
Omega::instance().scanPlugins(string(PREFIX "/lib/yade" SUFFIX "/gui" ));
+ Omega::instance().scanPlugins(string(PREFIX "/lib/yade" SUFFIX "/extra" ));
Omega::instance().init();
// make directory for temporaries
Added: trunk/extra/PeriodicInsertionSortCollider.cpp
===================================================================
--- trunk/extra/PeriodicInsertionSortCollider.cpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/extra/PeriodicInsertionSortCollider.cpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -0,0 +1,247 @@
+// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
+
+#include"PeriodicInsertionSortCollider.hpp"
+#include<yade/pkg-common/InsertionSortCollider.hpp>
+#include<yade/core/MetaBody.hpp>
+#include<yade/core/Interaction.hpp>
+#include<yade/core/InteractionContainer.hpp>
+#include<yade/pkg-common/BoundingVolumeMetaEngine.hpp>
+#include<yade/pkg-common/VelocityBins.hpp>
+#include<yade/pkg-dem/NewtonsDampedLaw.hpp>
+
+#include<algorithm>
+#include<vector>
+#include<boost/static_assert.hpp>
+
+using namespace std;
+
+YADE_PLUGIN((PeriodicInsertionSortCollider))
+CREATE_LOGGER(PeriodicInsertionSortCollider);
+
+Real PeriodicInsertionSortCollider::cellWrap(const Real x, const Real x0, const Real x1, long& period){
+ Real xNorm=(x-x0)/(x1-x0);
+ period=(long)floor(xNorm); // FIXME: some people say this is very slow
+ return x0+(xNorm-period)*(x1-x0);
+}
+
+
+// return true if bodies bb overlap in all 3 dimensions
+bool PeriodicInsertionSortCollider::spatialOverlap(body_id_t id1, body_id_t id2,MetaBody* rb, Vector3<int>& periods) const {
+ assert(id1!=id2) // programming error, or weird bodies (too large?)
+ for(int axis=0; axis<3; axis++){
+ Real dim=rb->cellMax[axis]-rb->cellMin[axis];
+ // wrap all 4 numbers to the period starting and the most minimal number
+ #if 0
+ Real mn=min(minima[3*id1+axis],minima[3*id2+axis])-0.001*dim; // avoid rounding issues
+ Real mx=max(maxima[3*id1+axis],maxima[3*id2+axis]);
+ TRVAR2(mn,mx);
+ #endif
+ // too big bodies in interaction
+ assert(maxima[3*id1+axis]-minima[3*id1+axis]<.99*dim); assert(maxima[3*id2+axis]-minima[3*id2+axis]<.99*dim);
+ // different way: find body of which when taken as period start will make the gap smaller
+ long p;
+ Real mn1w=cellWrap(minima[3*id1+axis],0,dim,p), mn2w=cellWrap(minima[3*id2+axis],0,dim,p);
+ Real wMn=(abs(mn2w-mn1w)<dim/2 ? mn1w : mn2w) -/*avoid rounding issues*/1e-4*dim; /* selected wrap base */
+ //TRVAR3(id1,id2,axis);
+ //TRVAR4(minima[3*id1+axis],maxima[3*id1+axis],minima[3*id2+axis],maxima[3*id2+axis]);
+ //TRVAR3(mn1w,mn2w,wMn);
+ long pmn1,pmx1,pmn2,pmx2;
+ Real mn1=cellWrap(minima[3*id1+axis],wMn,wMn+dim,pmn1), mx1=cellWrap(maxima[3*id1+axis],wMn,wMn+dim,pmx1);
+ Real mn2=cellWrap(minima[3*id2+axis],wMn,wMn+dim,pmn2), mx2=cellWrap(maxima[3*id2+axis],wMn,wMn+dim,pmx2);
+ //TRVAR4(mn1,mx1,mn2,mx2);
+ //TRVAR4(pmn1,pmx1,pmn2,pmx2);
+ assert(pmn1==pmx1); assert(pmn2==pmx2);
+ periods[axis]=(int)(pmn1-pmn2);
+ if(!(mn1<=mx2 && mx1 >= mn2)) return false;
+ }
+ //LOG_TRACE("Returning true for #"<<id1<<"+#"<<id2<<", periods "<<periods);
+ return true;
+}
+
+// called by the insertion sort if 2 bodies swapped their bounds
+void PeriodicInsertionSortCollider::handleBoundInversion(body_id_t id1, body_id_t id2, InteractionContainer* interactions, MetaBody* rb){
+ // do bboxes overlap in all 3 dimensions?
+ Vector3<int> periods;
+ bool overlap=spatialOverlap(id1,id2,rb,periods);
+ // existing interaction?
+ const shared_ptr<Interaction>& I=interactions->find(id1,id2);
+ bool hasInter=(bool)I;
+ // interaction doesn't exist and shouldn't, or it exists and should
+ if(!overlap && !hasInter) return;
+ if(overlap && hasInter){ return; }
+ // create interaction if not yet existing
+ if(overlap && !hasInter){ // second condition only for readability
+ if(!Collider::mayCollide(Body::byId(id1,rb).get(),Body::byId(id2,rb).get())) return;
+ // LOG_TRACE("Creating new interaction #"<<id1<<"+#"<<id2);
+ shared_ptr<Interaction> newI=shared_ptr<Interaction>(new Interaction(id1,id2));
+ newI->cellDist=periods;
+ interactions->insert(newI);
+ return;
+ }
+ if(!overlap && hasInter){ if(!I->isReal()) interactions->erase(id1,id2); return; }
+ assert(false); // unreachable
+}
+
+void PeriodicInsertionSortCollider::insertionSort(VecBounds& v, InteractionContainer* interactions, MetaBody* rb, bool doCollide){
+ // It seems that due to wrapping, it is not predetermined, how many times should we traverse the container
+ // We therefore add one blank traversal at the end that finds no inversions; then stop
+ long &loIdx=v.loIdx, &size=v.size; // shorthands
+ for(bool hadInversion=true; hadInversion; ){
+ hadInversion=false;
+ long cnt=0; //loIdx
+ for(long i=v.norm(loIdx-1); cnt++<size; i=v.norm(i-1)){
+ long i_1=v.norm(i-1), loIdx_1=v.norm(loIdx-1);
+ // fast test, if the first pair is inverted
+ if(v[i].coord<v[i_1].coord-(i_1==loIdx_1 ? v.cellDim : 0) ){
+ // v.dump(cerr);
+ hadInversion=true; Bound vi=v[i]; int j; const bool viBB=vi.flags.hasBB;
+ for(j=i_1; vi.coord<v[j].coord-(j==v.norm(loIdx-1) ? v.cellDim : 0); j=v.norm(j-1)) {
+ //{ Bound vj1=v[v.norm(j+1)]; v[v.norm(j+1)]=vi;
+ //v[v.norm(j+1)]=vj1; }
+ long j1=v.norm(j+1); // j2=v.norm(j+2);
+ //LOG_TRACE("Inversion of i="<<i<<"(#"<<vi.id<<" @ "<<vi.coord<<") j="<<j<<"(#"<<v[j].id<<" @ "<<v[j].coord<<"); j1="<<j1<<", j2="<<j2); v.dump(cerr);
+ v[j1]=v[j];
+ //if(v[j1].coord>v.cellMax && j2==loIdx){ v[j1].period+=1; v[j1].coord-=v.cellDim; loIdx=v.norm(loIdx-1); }
+ if(j1==loIdx) { assert(v[j1].coord>v.cellMax); v[j1].period+=1; v[j1].coord-=v.cellDim; loIdx=v.norm(loIdx-1); }
+ else if (vi.coord<v.cellMin && j==loIdx){ vi.period-=1; vi.coord+=v.cellDim; loIdx=v.norm(loIdx+1); }
+ if(doCollide && viBB && v[j].flags.hasBB) handleBoundInversion(vi.id,v[j].id,interactions,rb);
+ //v.dump(cerr);
+ }
+ v[v.norm(j+1)]=vi;
+ //LOG_TRACE("DONE:"); v.dump(cerr);
+ }
+ }
+ }
+}
+
+void PeriodicInsertionSortCollider::VecBounds::update(MetaBody* rb, int axis){
+ assert(axis>=0 && axis<3);
+ if(!rb->isPeriodic){
+ // maybe just set cell from -inf to +inf and go ahead?
+ LOG_FATAL("MetaBody::isPeriodic is false, unable to continue!"); throw;
+ }
+ cellMin=rb->cellMin[axis]; cellMax=rb->cellMax[axis]; cellDim=cellMax-cellMin;
+ size=vec.size();
+}
+
+void PeriodicInsertionSortCollider::action(MetaBody* rb){
+ long nBodies=(long)rb->bodies->size();
+ InteractionContainer* interactions=rb->interactions.get();
+
+ // pre-conditions
+ // adjust storage size
+ bool doInitSort=false;
+ if(XX.size!=2*nBodies){
+ LOG_DEBUG("Resize bounds containers from "<<XX.size<<" to "<<nBodies*2<<", will std::sort.");
+ // bodies deleted; clear the container completely, and do as if all bodies were added (rather slow…)
+ // future possibility: insertion sort with such operator that deleted bodies would all go to the end, then just trim bounds
+ if(2*nBodies<XX.size){ XX.vec.clear(); YY.vec.clear(); ZZ.vec.clear(); }
+ // more than 100 bodies was added, do initial sort again
+ // maybe: should rather depend on ratio of added bodies to those already present...?
+ if(2*nBodies-XX.size>200 || XX.size==0) doInitSort=true;
+ XX.vec.reserve(2*nBodies); YY.vec.reserve(2*nBodies); ZZ.vec.reserve(2*nBodies);
+ assert((XX.vec.size()%2)==0);
+ for(size_t id=XX.vec.size()/2; id<(size_t)nBodies; id++){
+ // add lower and upper bounds; coord is not important, will be updated from bb shortly
+ XX.vec.push_back(Bound(0,id,/*isMin=*/true)); XX.vec.push_back(Bound(0,id,/*isMin=*/false));
+ YY.vec.push_back(Bound(0,id, true)); YY.vec.push_back(Bound(0,id, false));
+ ZZ.vec.push_back(Bound(0,id, true)); ZZ.vec.push_back(Bound(0,id, false));
+ }
+ // XX.dump(cerr); YY.dump(cerr); ZZ.dump(cerr);
+ }
+ if(minima.size()!=3*(size_t)nBodies){ minima.resize(3*nBodies); maxima.resize(3*nBodies); }
+ assert(XX.vec.size()==2*rb->bodies->size());
+ //PERI: copy periodiciy information
+ XX.update(rb,0); YY.update(rb,1); ZZ.update(rb,2);
+
+ // copy bounds along given axis into our arrays
+ for(size_t i=0; i<2*(size_t)nBodies; i++){
+ const body_id_t& idXX=XX[i].id; const body_id_t& idYY=YY[i].id; const body_id_t& idZZ=ZZ[i].id;
+ const shared_ptr<BoundingVolume>& bvXX=Body::byId(idXX,rb)->boundingVolume; const shared_ptr<BoundingVolume>& bvYY=Body::byId(idYY,rb)->boundingVolume; const shared_ptr<BoundingVolume>& bvZZ=Body::byId(idZZ,rb)->boundingVolume;
+ // copy bounds from boundingVolume if there is one, otherwise use position; store what was used in the flags.hasBB bit
+ // PERI: add current period number to the coordinate
+ XX[i].coord=((XX[i].flags.hasBB=(bool)bvXX) ? (XX[i].flags.isMin ? bvXX->min[0] : bvXX->max[0]) : (Body::byId(idXX,rb)->physicalParameters->se3.position[0])) - XX.cellDim*XX[i].period;
+ YY[i].coord=((YY[i].flags.hasBB=(bool)bvYY) ? (YY[i].flags.isMin ? bvYY->min[1] : bvYY->max[1]) : (Body::byId(idYY,rb)->physicalParameters->se3.position[1])) - YY.cellDim*YY[i].period;
+ ZZ[i].coord=((ZZ[i].flags.hasBB=(bool)bvZZ) ? (ZZ[i].flags.isMin ? bvZZ->min[2] : bvZZ->max[2]) : (Body::byId(idZZ,rb)->physicalParameters->se3.position[2])) - ZZ.cellDim*ZZ[i].period;
+ // and for each body, copy its minima and maxima arrays as well
+ if(XX[i].flags.isMin){
+ BOOST_STATIC_ASSERT(sizeof(Vector3r)==3*sizeof(Real));
+ if(bvXX) {
+ memcpy(&minima[3*idXX],&bvXX->min,3*sizeof(Real)); memcpy(&maxima[3*idXX],&bvXX->max,3*sizeof(Real)); // ⇐ faster than 6 assignments
+ }
+ else{ const Vector3r& pos=Body::byId(idXX,rb)->physicalParameters->se3.position; memcpy(&minima[3*idXX],pos,3*sizeof(Real)); memcpy(&maxima[3*idXX],pos,3*sizeof(Real)); }
+ // PERI: add periods, but such that both minimum and maximum is within the cell!
+ Vector3r period(XX[i].period*XX.cellDim,YY[i].period*YY.cellDim,ZZ[i].period*ZZ.cellDim);
+ *(Vector3r*)(&minima[3*idXX])+=period; *(Vector3r*)(&maxima[3*idXX])+=period; //ugh
+ }
+ }
+
+ // process interactions that the constitutive law asked to be erased
+ interactions->erasePending(*this,rb);
+ LOG_DEBUG("Step "<<Omega::instance().getCurrentIteration());
+ ZZ.dump(cerr);
+ // XX.dump(cerr); YY.dump(cerr); ZZ.dump(cerr);
+
+ // sort
+ if(!doInitSort){
+ /* each inversion in insertionSort calls handleBoundInversion, which in turns may add/remove interaction */
+ insertionSort(XX,interactions,rb); insertionSort(YY,interactions,rb); insertionSort(ZZ,interactions,rb);
+ }
+ else {
+ if(doInitSort){
+ // the initial sort is in independent in 3 dimensions, may be run in parallel
+ // it seems that there is no time gain running this in parallel, though
+ std::sort(XX.vec.begin(),XX.vec.end()); std::sort(YY.vec.begin(),YY.vec.end()); std::sort(ZZ.vec.begin(),ZZ.vec.end());
+ } else { // sortThenCollide
+ insertionSort(XX,interactions,rb,false); insertionSort(YY,interactions,rb,false); insertionSort(ZZ,interactions,rb,false);
+ }
+ // traverse the container along requested axis
+ assert(sortAxis==0 || sortAxis==1 || sortAxis==2);
+ VecBounds& V=(sortAxis==0?XX:(sortAxis==1?YY:ZZ));
+ // go through potential aabb collisions, create interactions as necessary
+ for(size_t i=0; i<2*(size_t)nBodies; i++){
+ // start from the lower bound
+ // skip bodies without bbox since we would possibly never meet the upper bound again (std::sort may not be stable) and we don't want to collide those anyway
+ if(!(V[i].flags.isMin && V[i].flags.hasBB)) continue;
+ const body_id_t& iid=V[i].id;
+ /* If std::sort swaps equal min/max bounds, there are 2 cases distinct cases to handle:
+
+ 1. i is inside V, j gets to the end of V (since it loops till the maxbound is found)
+ here, we swap min/max to get the in the right order and continue with the loop over i;
+ next time this j-bound is handled (with a different i, for sure), it will be OK.
+ 2. i is at the end of V, therefore (j=i+1)==2*nBodies, therefore V[j] doesn't exist (past the end)
+ here, we can just check for that and break the loop if it happens.
+ It is the last i that we process, nothing will come after.
+
+ NOTE: XX,YY,ZZ containers don't guarantee that i_min<i_max. This is needed only here and is
+ handled only for the sortAxis. Functionality-wise, this has no impact on further collision
+ detection, though.
+ */
+ // TRVAR3(i,iid,V[i].coord);
+ // go up until we meet the upper bound
+ for(size_t j=i+1; V[j].id!=iid && /* handle case 2. of swapped min/max */ j<2*(size_t)nBodies; j++){
+ const body_id_t& jid=V[j].id;
+ /// Not sure why this doesn't work. If this condition is commented out, we have exact same interactions as from SpatialQuickSort. Otherwise some interactions are missing!
+ // skip bodies with smaller (arbitrary, could be greater as well) id, since they will detect us when their turn comes
+ //if(jid<iid) { /* LOG_TRACE("Skip #"<<V[j].id<<(V[j].flags.isMin?"(min)":"(max)")<<" with "<<iid<<" (smaller id)"); */ continue; }
+ // take 2 of the same condition (only handle collision [min_i..max_i]+min_j, not [min_i..max_i]+min_i (symmetric)
+ if(!V[j].flags.isMin) continue;
+ /* abuse the same function here; since it does spatial overlap check first, it is OK to use it */
+ handleBoundInversion(iid,jid,interactions,rb);
+ // now we are at the last element, but we still have not met the upper bound of V[i].id
+ // that means that the upper bound is before the upper one; that can only happen if they
+ // are equal and the unstable std::sort has swapped them. In that case, we need to go reverse
+ // from V[i] until we meet the upper bound and swap the isMin flag
+ if(j==2*(size_t)nBodies-1){ /* handle case 1. of swapped min/max */
+ size_t k=i-1;
+ while(V[k].id!=iid && k>0) k--;
+ assert(V[k].id==iid); // if this fails, we didn't meet the other bound in the downwards sense either; that should never happen
+ assert(!V[k].flags.isMin); // the lower bound should be maximum in this (exceptional) case; we will fix that now
+ V[k].flags.isMin=true; V[i].flags.isMin=false;
+ LOG_DEBUG("Swapping coincident min/max of #"<<iid<<" at positions "<<k<<" and "<<i<<" (coords: "<<V[k].coord<<" and "<<V[i].coord<<")");
+ break; // would happen anyways
+ }
+ }
+ }
+ }
+}
Added: trunk/extra/PeriodicInsertionSortCollider.hpp
===================================================================
--- trunk/extra/PeriodicInsertionSortCollider.hpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/extra/PeriodicInsertionSortCollider.hpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -0,0 +1,60 @@
+// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
+
+#pragma once
+#include<yade/core/Collider.hpp>
+class InteractionContainer;
+
+class PeriodicInsertionSortCollider: public Collider{
+ //! struct for storing bounds of bodies
+ struct Bound{
+ //! coordinate along the given sortAxis
+ Real coord;
+ //! id of the body this bound belongs to
+ body_id_t id;
+ //! periodic cell coordinate
+ int period;
+ //! is it the minimum (true) or maximum (false) bound?
+ struct{ unsigned hasBB:1; unsigned isMin:1; } flags;
+ Bound(Real coord_, body_id_t id_, bool isMin): coord(coord_), id(id_), period(0){ flags.isMin=isMin; }
+ bool operator<(const Bound& b) const {return coord<b.coord;}
+ bool operator>(const Bound& b) const {return coord>b.coord;}
+ };
+ struct VecBounds{
+ std::vector<Bound> vec;
+ Real cellMin, cellMax, cellDim;
+ long size;
+ // index of the lowest coordinate element, before which the container wraps
+ long loIdx;
+ Bound& operator[](long idx){ assert(idx<size && idx>=0); return vec[idx]; }
+ const Bound& operator[](long idx) const { assert(idx<size && idx>=0); return vec[idx]; }
+ // update periodic properties and size from MetaBody
+ void update(MetaBody* rb, int axis);
+ // normalize given index to the right range (wraps around)
+ long norm(long i) const { if(i<0) i+=size; long ret=i%size; assert(ret>=0 && ret<size); return ret;}
+ VecBounds(): size(0), loIdx(0){}
+ void dump(ostream& os){ string ret; for(size_t i=0; i<vec.size(); i++) os<<((long)i==loIdx?"@@ ":"")<<vec[i].coord<<"(id="<<vec[i].id<<","<<(vec[i].flags.isMin?"min":"max")<<",p"<<vec[i].period<<") "; os<<endl;}
+ };
+ private:
+ //! storage for bounds
+ VecBounds XX,YY,ZZ;
+ //! storage for bb maxima and minima
+ std::vector<Real> maxima, minima;
+
+ void insertionSort(VecBounds& v,InteractionContainer*,MetaBody*,bool doCollide=true);
+ void handleBoundInversion(body_id_t,body_id_t,InteractionContainer*,MetaBody*);
+ bool spatialOverlap(body_id_t,body_id_t, MetaBody*, Vector3<int>&) const;
+ static Real cellWrap(Real,Real,Real,long&);
+
+ public:
+ //! axis for the initial sort
+ int sortAxis;
+ //! Predicate called from loop within InteractionContainer::erasePending
+ bool shouldBeErased(body_id_t id1, body_id_t id2, MetaBody* rb) const { Vector3<int> foo; return !spatialOverlap(id1,id2,rb,foo); }
+
+ PeriodicInsertionSortCollider(): sortAxis(0) { }
+ virtual void action(MetaBody*);
+ REGISTER_CLASS_AND_BASE(PeriodicInsertionSortCollider,Collider);
+ REGISTER_ATTRIBUTES(Collider,(sortAxis));
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(PeriodicInsertionSortCollider);
Modified: trunk/extra/SConscript
===================================================================
--- trunk/extra/SConscript 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/extra/SConscript 2009-08-06 22:05:28 UTC (rev 1928)
@@ -1,10 +1,12 @@
# vim: set filetype=python :
Import('*')
+linkPlugins=env['linkPlugins']
+
import os.path, os
env.Install('$PREFIX/lib/yade$SUFFIX/extra',[
-
+ env.SharedLibrary('PeriodicInsertionSortCollider',['PeriodicInsertionSortCollider.cpp'],LIBS=env['LIBS']+linkPlugins(['InsertionSortCollider'])),
])
Modified: trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -76,7 +76,7 @@
if(XX.size()!=2*rb->bodies->size()) return true;
// we wouldn't run in this step; in that case, just delete pending interactions
// this is done in ::action normally, but it would make the call counters not reflect the stride
- rb->interactions->erasePending(*this);
+ rb->interactions->erasePending(*this,rb);
return false;
}
#endif
@@ -182,7 +182,7 @@
ISC_CHECKPOINT("copy");
// process interactions that the constitutive law asked to be erased
- interactions->erasePending(*this);
+ interactions->erasePending(*this,rb);
ISC_CHECKPOINT("erase");
// sort
Modified: trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -86,7 +86,7 @@
// This makes the collider non-persistent, not remembering last state
bool sortThenCollide;
//! Predicate called from loop within InteractionContainer::erasePending
- bool shouldBeErased(body_id_t id1, body_id_t id2) const { return !spatialOverlap(id1,id2); }
+ bool shouldBeErased(body_id_t id1, body_id_t id2, MetaBody*) const { return !spatialOverlap(id1,id2); }
#ifdef COLLIDE_STRIDED
virtual bool isActivated(MetaBody*);
#endif
Modified: trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -84,7 +84,7 @@
// timingDeltas->checkpoint("minMaxUpdate");
- ncb->interactions->erasePending(*this);
+ ncb->interactions->erasePending(*this,ncb);
// timingDeltas->checkpoint("deleteInvalid");
Modified: trunk/py/yadeWrapper/yadeWrapper.cpp
===================================================================
--- trunk/py/yadeWrapper/yadeWrapper.cpp 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/py/yadeWrapper/yadeWrapper.cpp 2009-08-06 22:05:28 UTC (rev 1928)
@@ -446,6 +446,16 @@
oa << YadeSimulation;
}
#endif
+ void periodicCell_set(python::tuple& t){
+ if(python::len(t)==2){ OMEGA.getRootBody()->cellMin=python::extract<Vector3r>(t[0]); OMEGA.getRootBody()->cellMax=python::extract<Vector3r>(t[1]); OMEGA.getRootBody()->isPeriodic=true; }
+ else if (python::len(t)==0) {OMEGA.getRootBody()->isPeriodic=false; }
+ else throw invalid_argument("Must pass either 2 Vector3's to set periodic cell, or () to disable periodicity (got "+lexical_cast<string>(python::len(t))+" instead).");
+ }
+ python::tuple periodicCell_get(){
+ vector<Vector3r> ret;
+ if(OMEGA.getRootBody()->isPeriodic){ return python::make_tuple(OMEGA.getRootBody()->cellMin,OMEGA.getRootBody()->cellMax); }
+ return python::make_tuple();
+ }
};
BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(omega_run_overloads,run,0,2);
BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(omega_saveTmp_overloads,saveTmp,0,1);
@@ -562,6 +572,7 @@
long Interaction_getId1(const shared_ptr<Interaction>& i){ return (long)i->getId1(); }
long Interaction_getId2(const shared_ptr<Interaction>& i){ return (long)i->getId2(); }
+python::tuple Interaction_getCellDist(const shared_ptr<Interaction>& i){ return python::make_tuple(i->cellDist[0],i->cellDist[1],i->cellDist[2]); }
void FileGenerator_generate(const shared_ptr<FileGenerator>& fg, string outFile){ fg->setFileName(outFile); fg->setSerializationLibrary("XMLFormatManager"); bool ret=fg->generateAndSave(); LOG_INFO((ret?"SUCCESS:\n":"FAILURE:\n")<<fg->message); if(ret==false) throw runtime_error("Generator reported error: "+fg->message); };
void FileGenerator_load(const shared_ptr<FileGenerator>& fg){ string xml(Omega::instance().tmpFilename()+".xml.bz2"); LOG_DEBUG("Using temp file "<<xml); FileGenerator_generate(fg,xml); pyOmega().load(xml); }
@@ -625,6 +636,7 @@
.add_property("timingEnabled",&pyOmega::timingEnabled_get,&pyOmega::timingEnabled_set,"Globally enable/disable timing services (see documentation of yade.timing).")
.add_property("bexSyncCount",&pyOmega::bexSyncCount_get,&pyOmega::bexSyncCount_set,"Counter for number of syncs in BexContainer, for profiling purposes.")
.add_property("numThreads",&pyOmega::numThreads_get /* ,&pyOmega::numThreads_set*/ ,"Get maximum number of threads openMP can use.")
+ .add_property("periodicCell",&pyOmega::periodicCell_get,&pyOmega::periodicCell_set, "Get/set periodic cell minimum and maximum (tuple of 2 Vector3's), or () for no periodicity.");
#ifdef YADE_BOOST_SERIALIZATION
.def("saveXML",&pyOmega::saveXML,"[EXPERIMENTAL] function saving to XML file using boost::serialization.")
#endif
@@ -761,7 +773,8 @@
.def_readwrite("geom",&Interaction::interactionGeometry)
.add_property("id1",&Interaction_getId1)
.add_property("id2",&Interaction_getId2)
- .add_property("isReal",&Interaction::isReal);
+ .add_property("isReal",&Interaction::isReal)
+ .add_property("cellDist",&Interaction_getCellDist);
EXPOSE_CXX_CLASS(InteractionPhysics);
EXPOSE_CXX_CLASS(InteractionGeometry);
EXPOSE_CXX_CLASS(FileGenerator)
Added: trunk/scripts/test/periodic-simple.py
===================================================================
--- trunk/scripts/test/periodic-simple.py 2009-08-06 08:40:14 UTC (rev 1927)
+++ trunk/scripts/test/periodic-simple.py 2009-08-06 22:05:28 UTC (rev 1928)
@@ -0,0 +1,24 @@
+from yade import log
+log.setLevel("PeriodicInsertionSortCollider",log.TRACE)
+O.engines=[
+ BexResetter(),
+ BoundingVolumeMetaEngine([InteractingSphere2AABB(),MetaInteractingGeometry2AABB()]),
+ PeriodicInsertionSortCollider(label='collider'),
+ InteractionDispatchers(
+ [ef2_Sphere_Sphere_Dem3DofGeom()],
+ [SimpleElasticRelationships()],
+ [Law2_Dem3Dof_Elastic_Elastic()],
+ ),
+ GravityEngine(gravity=[0,0,-10]),
+ NewtonsDampedLaw(damping=.1)
+]
+O.bodies.append(utils.sphere([0,0,2],1,dynamic=False,density=1000))
+O.bodies.append(utils.sphere([0,0,3],1,density=1000))
+O.periodicCell=((-5,-5,0),(5,5,10))
+O.dt=utils.PWaveTimeStep()
+O.saveTmp()
+from yade import qt
+qt.Controller()
+qt.View()
+O.run(17)
+
