yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01624
[svn] r1929 - in trunk: extra gui/qt3 pkg/common/Engine/MetaEngine pkg/common/Engine/StandAloneEngine pkg/common/RenderingEngine/OpenGLRenderingEngine scripts/test
Author: eudoxos
Date: 2009-08-07 12:27:49 +0200 (Fri, 07 Aug 2009)
New Revision: 1929
Modified:
trunk/extra/PeriodicInsertionSortCollider.cpp
trunk/extra/PeriodicInsertionSortCollider.hpp
trunk/gui/qt3/GLViewer.cpp
trunk/gui/qt3/GLViewer.hpp
trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp
trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp
trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
trunk/pkg/common/RenderingEngine/OpenGLRenderingEngine/OpenGLRenderingEngine.cpp
trunk/pkg/common/RenderingEngine/OpenGLRenderingEngine/OpenGLRenderingEngine.hpp
trunk/scripts/test/periodic-simple.py
Log:
1. Beta version of periodic boundary conditions (try scripts/test/periodic-simple.py
Modified: trunk/extra/PeriodicInsertionSortCollider.cpp
===================================================================
--- trunk/extra/PeriodicInsertionSortCollider.cpp 2009-08-06 22:05:28 UTC (rev 1928)
+++ trunk/extra/PeriodicInsertionSortCollider.cpp 2009-08-07 10:27:49 UTC (rev 1929)
@@ -18,34 +18,33 @@
YADE_PLUGIN((PeriodicInsertionSortCollider))
CREATE_LOGGER(PeriodicInsertionSortCollider);
-Real PeriodicInsertionSortCollider::cellWrap(const Real x, const Real x0, const Real x1, long& period){
+Real PeriodicInsertionSortCollider::cellWrap(const Real x, const Real x0, const Real x1, int& period){
Real xNorm=(x-x0)/(x1-x0);
- period=(long)floor(xNorm); // FIXME: some people say this is very slow
+ period=(int)floor(xNorm); // FIXME: some people say this is very slow
return x0+(xNorm-period)*(x1-x0);
}
+Real PeriodicInsertionSortCollider::cellWrapRel(const Real x, const Real x0, const Real x1){
+ Real xNorm=(x-x0)/(x1-x0);
+ return (xNorm-floor(xNorm))*(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?)
+ 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 */
+ Real m1=minima[3*id1+axis],m2=minima[3*id2+axis];
+ Real wMn=(cellWrapRel(m1,m2,m2+dim)<cellWrapRel(m2,m1,m1+dim)) ? m2 : m1;
//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;
+ //TRVAR2(cellWrapRel(m1,m2,m2+dim),cellWrapRel(m2,m1,m1+dim));
+ //TRVAR3(m1,m2,wMn);
+ int 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);
@@ -93,16 +92,17 @@
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);
+ //v.dump(cerr);
+ if(i==loIdx && v[i].coord<v.cellMin){ v[i].period-=1; v[i].coord+=v.cellDim; loIdx=v.norm(loIdx+1); }
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);
+ long j1=v.norm(j+1);
+ //LOG_TRACE("Inversion of i="<<i<<"(#"<<vi.id<<" @ "<<vi.coord<<") j="<<j<<"(#"<<v[j].id<<" @ "<<v[j].coord<<"); j1="<<j1); 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); }
+ 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);
@@ -163,6 +163,13 @@
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;
+ // PERI: at the initial step, fix periods of bodies
+ // doInitSort is also called when bodies are just added; changing the period should not have influence here, though.
+ if(doInitSort){
+ if(XX[i].coord<XX.cellMin || XX[i].coord>=XX.cellMax) XX[i].coord=cellWrap(XX[i].coord,XX.cellMin,XX.cellMax,XX[i].period);
+ if(YY[i].coord<XX.cellMin || YY[i].coord>=YY.cellMax) YY[i].coord=cellWrap(YY[i].coord,YY.cellMin,YY.cellMax,YY[i].period);
+ if(ZZ[i].coord<ZZ.cellMin || ZZ[i].coord>=ZZ.cellMax) ZZ[i].coord=cellWrap(ZZ[i].coord,ZZ.cellMin,ZZ.cellMax,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));
@@ -171,15 +178,15 @@
}
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
+ //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);
+ //LOG_TRACE("Step "<<Omega::instance().getCurrentIteration());
+ //ZZ.dump(cerr);
// XX.dump(cerr); YY.dump(cerr); ZZ.dump(cerr);
// sort
@@ -188,13 +195,10 @@
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);
- }
+ // 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());
+
// traverse the container along requested axis
assert(sortAxis==0 || sortAxis==1 || sortAxis==2);
VecBounds& V=(sortAxis==0?XX:(sortAxis==1?YY:ZZ));
@@ -219,7 +223,7 @@
*/
// 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++){
+ for(size_t j=i+1; /* handle case 2. of swapped min/max */ j<2*(size_t)nBodies && V[j].id!=iid; 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
Modified: trunk/extra/PeriodicInsertionSortCollider.hpp
===================================================================
--- trunk/extra/PeriodicInsertionSortCollider.hpp 2009-08-06 22:05:28 UTC (rev 1928)
+++ trunk/extra/PeriodicInsertionSortCollider.hpp 2009-08-07 10:27:49 UTC (rev 1929)
@@ -43,8 +43,10 @@
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&);
+ static Real cellWrap(const Real, const Real, const Real, int&);
+ static Real cellWrapRel(const Real, const Real, const Real);
+
public:
//! axis for the initial sort
int sortAxis;
Modified: trunk/gui/qt3/GLViewer.cpp
===================================================================
--- trunk/gui/qt3/GLViewer.cpp 2009-08-06 22:05:28 UTC (rev 1928)
+++ trunk/gui/qt3/GLViewer.cpp 2009-08-07 10:27:49 UTC (rev 1929)
@@ -323,6 +323,18 @@
else if(e->key()!=Qt::Key_Escape && e->key()!=Qt::Key_Space) QGLViewer::keyPressEvent(e);
updateGL();
}
+/* Center the scene such that periodic cell is contained in the view */
+void GLViewer::centerPeriodic(){
+ MetaBody* rb=Omega::instance().getRootBody().get();
+ assert(rb->isPeriodic);
+ Vector3r center=.5*(rb->cellMin+rb->cellMax);
+ Vector3r halfSize=.5*(rb->cellMax-rb->cellMin);
+ float radius=std::max(halfSize[0],std::max(halfSize[1],halfSize[2]));
+ LOG_DEBUG("Periodic scene center="<<center<<", halfSize="<<halfSize<<", radius="<<radius);
+ setSceneCenter(qglviewer::Vec(center[0],center[1],center[2]));
+ setSceneRadius(radius*1.5);
+ showEntireScene();
+}
/* Calculate medians for x, y and z coordinates of all bodies;
*then set scene center to median position and scene radius to 2*inter-quartile distance.
@@ -332,6 +344,7 @@
* "central" (where most bodies is) part very small or even invisible.
*/
void GLViewer::centerMedianQuartile(){
+ if(Omega::instance().getRootBody()->isPeriodic){ centerPeriodic(); return; }
long nBodies=Omega::instance().getRootBody()->bodies->size();
if(nBodies<4) {
LOG_INFO("Less than 4 bodies, median makes no sense; calling centerScene() instead.");
@@ -357,6 +370,7 @@
void GLViewer::centerScene(){
MetaBody* rb=Omega::instance().getRootBody().get();
if (!rb) return;
+ if(rb->isPeriodic){ centerPeriodic(); return; }
if(rb->bodies->size()<renderer->selectBodyLimit){LOG_INFO("Less than "+lexical_cast<string>(renderer->selectBodyLimit)+" bodies, moving possible. Select with shift, press 'm' to move.");}
else{LOG_INFO("More than "+lexical_cast<string>(renderer->selectBodyLimit)+" (OpenGLRenderingEngine::selectBodyLimit) bodies. Moving not possible.");}
Modified: trunk/gui/qt3/GLViewer.hpp
===================================================================
--- trunk/gui/qt3/GLViewer.hpp 2009-08-06 22:05:28 UTC (rev 1928)
+++ trunk/gui/qt3/GLViewer.hpp 2009-08-07 10:27:49 UTC (rev 1929)
@@ -66,6 +66,7 @@
virtual void draw();
virtual void drawWithNames();
void centerScene();
+ void centerPeriodic();
void mouseMovesCamera();
void mouseMovesManipulatedFrame(qglviewer::Constraint* c=NULL);
void resetManipulation();
Modified: trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp 2009-08-06 22:05:28 UTC (rev 1928)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp 2009-08-07 10:27:49 UTC (rev 1929)
@@ -34,6 +34,7 @@
LOG_WARN("Interactions pending erase found (erased), no collider being used?");
alreadyWarnedNoCollider=true;
}
+ Vector3r cellSize; if(rootBody->isPeriodic) cellSize=rootBody->cellMax-rootBody->cellMin;
#ifdef YADE_OPENMP
const long size=rootBody->interactions->size();
#pragma omp parallel for schedule(guided)
@@ -72,7 +73,12 @@
assert(I->functorCache.geom);
bool wasReal=I->isReal();
- bool geomCreated=I->functorCache.geom->go(b1->interactingGeometry,b2->interactingGeometry,b1->physicalParameters->se3, b2->physicalParameters->se3,I);
+ bool geomCreated;
+ if(!rootBody->isPeriodic) geomCreated=I->functorCache.geom->go(b1->interactingGeometry,b2->interactingGeometry,b1->physicalParameters->se3, b2->physicalParameters->se3,I);
+ else{ // handle periodicity
+ Se3r se32=b2->physicalParameters->se3; se32.position+=Vector3r(I->cellDist[0]*cellSize[0],I->cellDist[1]*cellSize[1],I->cellDist[2]*cellSize[2]);
+ geomCreated=I->functorCache.geom->go(b1->interactingGeometry,b2->interactingGeometry,b1->physicalParameters->se3,se32,I);
+ }
if(!geomCreated){
if(wasReal) rootBody->interactions->requestErase(I->getId1(),I->getId2()); // fully created interaction without geometry is reset and perhaps erased in the next step
continue; // in any case don't care about this one anymore
Modified: trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp 2009-08-06 22:05:28 UTC (rev 1928)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp 2009-08-07 10:27:49 UTC (rev 1929)
@@ -51,6 +51,7 @@
}
shared_ptr<BodyContainer>& bodies = ncb->bodies;
+ Vector3r cellSize; if(ncb->isPeriodic) cellSize=ncb->cellMax-ncb->cellMin;
#ifdef YADE_OPENMP
const long size=ncb->transientInteractions->size();
#pragma omp parallel for
@@ -63,10 +64,16 @@
const shared_ptr<Body>& b2=(*bodies)[I->getId2()];
bool wasReal=I->isReal();
if (!b1->interactingGeometry || !b2->interactingGeometry) { assert(!wasReal); continue; } // some bodies do not have interactingGeometry
- bool geomCreated=operator()(b1->interactingGeometry, b2->interactingGeometry, b1->physicalParameters->se3, b2->physicalParameters->se3, I);
+ bool geomCreated;
+ if(!ncb->isPeriodic){
+ geomCreated=operator()(b1->interactingGeometry, b2->interactingGeometry, b1->physicalParameters->se3, b2->physicalParameters->se3, I);
+ } else{
+ Se3r se32=b2->physicalParameters->se3; se32.position+=Vector3r(I->cellDist[0]*cellSize[0],I->cellDist[1]*cellSize[1],I->cellDist[2]*cellSize[2]); // add periodicity to the position of the 2nd body
+ geomCreated=operator()(b1->interactingGeometry, b2->interactingGeometry, b1->physicalParameters->se3, se32, I);
+ }
// reset && erase interaction that existed but now has no geometry anymore
if(wasReal && !geomCreated){ ncb->interactions->requestErase(I->getId1(),I->getId2()); }
}
}
-YADE_PLUGIN((InteractionGeometryMetaEngine));
\ No newline at end of file
+YADE_PLUGIN((InteractionGeometryMetaEngine));
Modified: trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp 2009-08-06 22:05:28 UTC (rev 1928)
+++ trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp 2009-08-07 10:27:49 UTC (rev 1929)
@@ -230,7 +230,7 @@
*/
// 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*nBodies; j++){
+ for(size_t j=i+1; /* handle case 2. of swapped min/max */ j<2*nBodies && V[j].id!=iid; 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
Modified: trunk/pkg/common/RenderingEngine/OpenGLRenderingEngine/OpenGLRenderingEngine.cpp
===================================================================
--- trunk/pkg/common/RenderingEngine/OpenGLRenderingEngine/OpenGLRenderingEngine.cpp 2009-08-06 22:05:28 UTC (rev 1928)
+++ trunk/pkg/common/RenderingEngine/OpenGLRenderingEngine/OpenGLRenderingEngine.cpp 2009-08-07 10:27:49 UTC (rev 1929)
@@ -109,16 +109,28 @@
numBodiesWhenRefSe3LastSet=rootBody->bodies->size();
numIterWhenRefSe3LastSet=Omega::instance().getCurrentIteration();
}
+/* mostly copied from PeriodicInsertionSortCollider
+ FIXME: common implementation somewhere */
+Real OpenGLRenderingEngine::wrapCell(const Real x, const Real x0, const Real x1){
+ Real xNorm=(x-x0)/(x1-x0);
+ return x0+(xNorm-floor(xNorm))*(x1-x0);
+}
+Vector3r OpenGLRenderingEngine::wrapCellPt(const Vector3r& pt, MetaBody* rb){
+ if(!rb->isPeriodic) return pt;
+ return Vector3r(wrapCell(pt[0],rb->cellMin[0],rb->cellMax[0]),wrapCell(pt[1],rb->cellMin[1],rb->cellMax[1]),wrapCell(pt[2],rb->cellMin[2],rb->cellMax[2]));
+}
+
void OpenGLRenderingEngine::setBodiesDispSe3(const shared_ptr<MetaBody>& rootBody){
FOREACH(const shared_ptr<Body>& b, *rootBody->bodies){
if(!b->physicalParameters) continue;
const Se3r& se3=b->physicalParameters->se3; const Se3r& refSe3=b->physicalParameters->refSe3; Se3r& dispSe3=b->physicalParameters->dispSe3;
- b->physicalParameters->isDisplayed=!pointClipped(se3.position);
+ Vector3r posCell=wrapCellPt(se3.position,rootBody.get());
+ b->physicalParameters->isDisplayed=!pointClipped(posCell);
// if no scaling, return quickly
- if(!(scaleDisplacements||scaleRotations)){ b->physicalParameters->dispSe3=b->physicalParameters->se3; continue; }
+ if(!(scaleDisplacements||scaleRotations||rootBody->isPeriodic)){ b->physicalParameters->dispSe3=b->physicalParameters->se3; continue; }
// apply scaling
- dispSe3.position=(scaleDisplacements ? diagMult(displacementScale,se3.position-refSe3.position)+refSe3.position : se3.position );
+ dispSe3.position=(scaleDisplacements ? diagMult(displacementScale,se3.position-refSe3.position)+wrapCellPt(refSe3.position,rootBody.get()) : posCell );
if(scaleRotations){
Quaternionr relRot=refSe3.orientation.Conjugate()*se3.orientation;
Vector3r axis; Real angle; relRot.ToAxisAngle(axis,angle);
@@ -127,7 +139,19 @@
} else {dispSe3.orientation=se3.orientation;}
}
}
+// draw periodic cell, if active
+void OpenGLRenderingEngine::drawPeriodicCell(MetaBody* rootBody){
+ if(!rootBody->isPeriodic) return;
+ glPushMatrix();
+ glColor3v(Vector3r(1,1,0));
+ Vector3r cent=.5*(rootBody->cellMin+rootBody->cellMax); Vector3r size=rootBody->cellMax-rootBody->cellMin;
+ glTranslate(cent[0],cent[1],cent[2]); glScale(size[0],size[1],size[2]);
+ glutWireCube(1);
+ glPopMatrix();
+}
+
+
void OpenGLRenderingEngine::render(const shared_ptr<MetaBody>& rootBody, body_id_t selection /*FIXME: not sure. maybe a list of selections, or maybe bodies themselves should remember if they are selected? */) {
assert(glutInitDone);
@@ -167,6 +191,8 @@
// set displayed Se3 of body (scaling) and isDisplayed (clipping)
setBodiesDispSe3(rootBody);
+ drawPeriodicCell(rootBody.get());
+
if (Show_DOF || Show_ID) renderDOF_ID(rootBody);
if (Body_geometrical_model){
if (Cast_shadows){
@@ -401,49 +427,78 @@
if(rootBody->geometricalModel) geometricalModelDispatcher(rootBody->geometricalModel,rootBody->physicalParameters,Body_wire);
}
-void OpenGLRenderingEngine::renderGeometricalModel(const shared_ptr<MetaBody>& rootBody){
+void OpenGLRenderingEngine::renderGeometricalModel(const shared_ptr<MetaBody>& rootBody){
const GLfloat ambientColorSelected[4]={10.0,0.0,0.0,1.0};
const GLfloat ambientColorUnselected[4]={0.5,0.5,0.5,1.0};
- if((rootBody->geometricalModel || Draw_inside) && Draw_inside) {
- FOREACH(const shared_ptr<Body> b, *rootBody->bodies){
- if(b->geometricalModel && ((b->getGroupMask() & Draw_mask) || b->getGroupMask()==0)){
- if(b->physicalParameters && !b->physicalParameters->isDisplayed) continue;
- const Se3r& se3=b->physicalParameters->dispSe3;
- glPushMatrix();
- Real angle; Vector3r axis; se3.orientation.ToAxisAngle(axis,angle);
- glTranslatef(se3.position[0],se3.position[1],se3.position[2]);
- glRotatef(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]);
- if(current_selection==b->getId() || b->geometricalModel->highlight){
- glLightModelfv(GL_LIGHT_MODEL_AMBIENT,ambientColorSelected);
- glColorMaterial(GL_FRONT_AND_BACK,GL_EMISSION);
- const Vector3r& h(current_selection==b->getId() ? highlightEmission0 : highlightEmission1);
- glColor4(h[0],h[1],h[2],.2);
- glColorMaterial(GL_FRONT_AND_BACK,GL_DIFFUSE);
+ if(rootBody->geometricalModel) geometricalModelDispatcher(rootBody->geometricalModel,rootBody->physicalParameters,Body_wire);
+ if(!Draw_inside) return;
+ FOREACH(const shared_ptr<Body> b, *rootBody->bodies){
+ if(!b->geometricalModel || (!((b->getGroupMask() & Draw_mask) || b->getGroupMask()==0))) continue;
+ if(b->physicalParameters && !b->physicalParameters->isDisplayed) continue;
+ const Se3r& se3=b->physicalParameters->dispSe3;
+ glPushMatrix();
+ Real angle; Vector3r axis; se3.orientation.ToAxisAngle(axis,angle);
+ glTranslatef(se3.position[0],se3.position[1],se3.position[2]);
+ glRotatef(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]);
+ if(current_selection==b->getId() || b->geometricalModel->highlight){
+ glLightModelfv(GL_LIGHT_MODEL_AMBIENT,ambientColorSelected);
+ glColorMaterial(GL_FRONT_AND_BACK,GL_EMISSION);
+ const Vector3r& h(current_selection==b->getId() ? highlightEmission0 : highlightEmission1);
+ glColor4(h[0],h[1],h[2],.2);
+ glColorMaterial(GL_FRONT_AND_BACK,GL_DIFFUSE);
- geometricalModelDispatcher(b->geometricalModel,b->physicalParameters,Body_wire);
+ geometricalModelDispatcher(b->geometricalModel,b->physicalParameters,Body_wire);
- glLightModelfv(GL_LIGHT_MODEL_AMBIENT,ambientColorUnselected);
- glColorMaterial(GL_FRONT_AND_BACK,GL_EMISSION);
- glColor3v(Vector3r::ZERO);
- glColorMaterial(GL_FRONT_AND_BACK,GL_DIFFUSE);
- } else {
- geometricalModelDispatcher(b->geometricalModel,b->physicalParameters,Body_wire);
+ glLightModelfv(GL_LIGHT_MODEL_AMBIENT,ambientColorUnselected);
+ glColorMaterial(GL_FRONT_AND_BACK,GL_EMISSION);
+ glColor3v(Vector3r::ZERO);
+ glColorMaterial(GL_FRONT_AND_BACK,GL_DIFFUSE);
+ } else {
+ geometricalModelDispatcher(b->geometricalModel,b->physicalParameters,Body_wire);
+ }
+ glPopMatrix();
+ if(current_selection==b->getId() || b->geometricalModel->highlight){
+ if(!b->boundingVolume || Body_wire || b->geometricalModel->wire) GLUtils::GLDrawInt(b->getId(),se3.position);
+ else {
+ // move the label towards the camera by the bounding box so that it is not hidden inside the body
+ const Vector3r& mn=b->boundingVolume->min; const Vector3r& mx=b->boundingVolume->max; const Vector3r& p=se3.position;
+ Vector3r ext(viewDirection[0]>0?p[0]-mn[0]:p[0]-mx[0],viewDirection[1]>0?p[1]-mn[1]:p[1]-mx[1],viewDirection[2]>0?p[2]-mn[2]:p[2]-mx[2]); // signed extents towards the camera
+ Vector3r dr=-1.01*(viewDirection.Dot(ext)*viewDirection);
+ GLUtils::GLDrawInt(b->getId(),se3.position+dr,Vector3r::ONE);
+ }
+ }
+ // if the body goes over the cell margin, draw it in all other positions with wire
+ if(b->boundingVolume && rootBody->isPeriodic){
+ const Vector3r& cellMin(rootBody->cellMin); const Vector3r& cellMax(rootBody->cellMax); Vector3r cellSize=cellMax-cellMin;
+ Vector3<int> bodyPer,minPer,maxPer;
+ for(int i=0; i<3; i++){
+ bodyPer[i]=(int)floor((b->physicalParameters->se3.position[i]-cellMin[i])/cellSize[i]);
+ minPer[i]=(int)floor((b->boundingVolume->min[i]-cellMin[i])/cellSize[i]);
+ maxPer[i]=(int)floor((b->boundingVolume->max[i]-cellMin[i])/cellSize[i]);
+ //assert(bodyPer[i]<=maxPer[i]); assert(bodyPer[i]>=minPer[i]);
+ }
+ /* m is bitmask from 3 couples (0…64=2^6) */
+ for(int m=0; m<64; m++){
+ // any mask containing 00 couple is invalid
+ if((!(m&1) && (!(m&2))) || (!(m&4) && (!(m&8))) || (!(m&16) && (!(m&32)))) continue;
+ Vector3r pt(se3.position);
+ bool isInside=false;
+ for(int j=0; j<3; j++){
+ if(m&(1<<(2*j))) {
+ if(m&(1<<(2*j+1))) { if(bodyPer[j]>=maxPer[j]) {isInside=true; break; } pt[j]-=cellSize[j]; }
+ else { if(bodyPer[j]<=minPer[j]){ isInside=true; break; } pt[j]+=cellSize[j]; }
+ }
}
+ if(isInside) continue;
+ if(pt==se3.position) continue; // shouldn't happen, but it happens :-(
+ glPushMatrix();
+ glTranslatev(pt);
+ glRotatef(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]);
+ geometricalModelDispatcher(b->geometricalModel,b->physicalParameters,/*Body_wire*/ true);
glPopMatrix();
- if(current_selection==b->getId() || b->geometricalModel->highlight){
- if(!b->boundingVolume || Body_wire || b->geometricalModel->wire) GLUtils::GLDrawInt(b->getId(),se3.position);
- else {
- // move the label towards the camera by the bounding box so that it is not hidden inside the body
- const Vector3r& mn=b->boundingVolume->min; const Vector3r& mx=b->boundingVolume->max; const Vector3r& p=se3.position;
- Vector3r ext(viewDirection[0]>0?p[0]-mn[0]:p[0]-mx[0],viewDirection[1]>0?p[1]-mn[1]:p[1]-mx[1],viewDirection[2]>0?p[2]-mn[2]:p[2]-mx[2]); // signed extents towards the camera
- Vector3r dr=-1.01*(viewDirection.Dot(ext)*viewDirection);
- GLUtils::GLDrawInt(b->getId(),se3.position+dr,Vector3r::ONE);
- }
- }
}
}
}
- if(rootBody->geometricalModel) geometricalModelDispatcher(rootBody->geometricalModel,rootBody->physicalParameters,Body_wire);
}
@@ -509,10 +564,10 @@
const Se3r& se3=b->physicalParameters->dispSe3;
if(b->interactingGeometry && ((b->getGroupMask()&Draw_mask) || b->getGroupMask()==0)){
glPushMatrix();
- Real angle; Vector3r axis; se3.orientation.ToAxisAngle(axis,angle);
- glTranslatef(se3.position[0],se3.position[1],se3.position[2]);
- glRotatef(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]);
- interactingGeometryDispatcher(b->interactingGeometry,b->physicalParameters,Body_wire);
+ Real angle; Vector3r axis; se3.orientation.ToAxisAngle(axis,angle);
+ glTranslatef(se3.position[0],se3.position[1],se3.position[2]);
+ glRotatef(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]);
+ interactingGeometryDispatcher(b->interactingGeometry,b->physicalParameters,Body_wire);
glPopMatrix();
}
}
Modified: trunk/pkg/common/RenderingEngine/OpenGLRenderingEngine/OpenGLRenderingEngine.hpp
===================================================================
--- trunk/pkg/common/RenderingEngine/OpenGLRenderingEngine/OpenGLRenderingEngine.hpp 2009-08-06 22:05:28 UTC (rev 1928)
+++ trunk/pkg/common/RenderingEngine/OpenGLRenderingEngine/OpenGLRenderingEngine.hpp 2009-08-07 10:27:49 UTC (rev 1929)
@@ -39,6 +39,13 @@
Real normSaw(Real t, Real period){ Real xi=(t-period*((int)(t/period)))/period; /* normalized value, (0-1〉 */ return (xi<.5?2*xi:2-2*xi); }
Real normSquare(Real t, Real period){ Real xi=(t-period*((int)(t/period)))/period; /* normalized value, (0-1〉 */ return (xi<.5?0:1); }
+ //! wrap number to interval x0…x1
+ Real wrapCell(const Real x, const Real x0, const Real x1);
+ //! wrap point to inside MetaBody's cell (identity if !MetaBody::isPeriodic)
+ Vector3r wrapCellPt(const Vector3r& pt, MetaBody* rb);
+ void drawPeriodicCell(MetaBody*);
+
+
private :
DynLibDispatcher< InteractionGeometry , GLDrawInteractionGeometryFunctor, void , TYPELIST_5(const shared_ptr<InteractionGeometry>&, const shared_ptr<Interaction>& , const shared_ptr<Body>&, const shared_ptr<Body>&, bool) > interactionGeometryDispatcher;
DynLibDispatcher< InteractionPhysics , GLDrawInteractionPhysicsFunctor, void , TYPELIST_5(const shared_ptr<InteractionPhysics>& , const shared_ptr<Interaction>&, const shared_ptr<Body>&, const shared_ptr<Body>&, bool) > interactionPhysicsDispatcher;
Modified: trunk/scripts/test/periodic-simple.py
===================================================================
--- trunk/scripts/test/periodic-simple.py 2009-08-06 22:05:28 UTC (rev 1928)
+++ trunk/scripts/test/periodic-simple.py 2009-08-07 10:27:49 UTC (rev 1929)
@@ -10,15 +10,15 @@
[Law2_Dem3Dof_Elastic_Elastic()],
),
GravityEngine(gravity=[0,0,-10]),
- NewtonsDampedLaw(damping=.1)
+ TranslationEngine(translationAxis=(1,0,0),velocity=10,subscribedBodies=[0]),
+ NewtonsDampedLaw(damping=.4)
]
-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.bodies.append(utils.sphere([-4,0,11],2,dynamic=False,density=1000))
+O.bodies.append(utils.sphere([0,-1,5.5],2,density=1000))
+O.bodies.append(utils.sphere([0,2,5.5],2,density=2000))
O.periodicCell=((-5,-5,0),(5,5,10))
-O.dt=utils.PWaveTimeStep()
+O.dt=.1*utils.PWaveTimeStep()
O.saveTmp()
from yade import qt
qt.Controller()
qt.View()
-O.run(17)
-
Follow ups
