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 <eudo...@arcig.cz> + +#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 <eudo...@arcig.cz> + +#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) + _______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp