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)
-
_______________________________________________
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
