yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02600
[Branch ~yade-dev/yade/trunk] Rev 1878: 1. Add Scene* Functor::scene, update it from all dispatchers at every step.
------------------------------------------------------------
revno: 1878
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Wed 2009-12-09 17:43:25 +0100
message:
1. Add Scene* Functor::scene, update it from all dispatchers at every step.
2. Add NewtonsDampedLaw::homotheticCellResize (0,1,2 for disabled, velocity or position update) (not tested at all)
3. Start removing Scene* ncb and Scene* rootBody etc from Engine::action as I stumble accross them. Please do the same, Engine::scene can be always used instead.
modified:
core/Dispatcher.hpp
core/Functor.hpp
pkg/common/Engine/Dispatcher/BoundDispatcher.cpp
pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp
pkg/common/Engine/Dispatcher/InteractionDispatchers.hpp
pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp
pkg/common/Engine/Dispatcher/InteractionPhysicsDispatcher.cpp
pkg/common/Engine/Dispatcher/LawDispatcher.cpp
pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp
pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp
pkg/dem/PreProcessor/TriaxialTest.cpp
py/timing.py
py/yadeWrapper/yadeWrapper.cpp
scripts/test-sphere-facet.py
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.
=== modified file 'core/Dispatcher.hpp'
--- core/Dispatcher.hpp 2009-12-01 14:56:39 +0000
+++ core/Dispatcher.hpp 2009-12-09 16:43:25 +0000
@@ -31,6 +31,9 @@
virtual void add(Functor*) {throw;}
virtual void add(string) {throw;}
+ //! Should be called by dispatcher before every loop (could be somehow optimized, since it will change very rarely)
+ void updateScenePtr() { FOREACH(shared_ptr<Functor> f, functorArguments){ f->scene=scene; } }
+
void storeFunctorName(const string& baseClassName1, const string& libName, shared_ptr<Functor> eu);
void storeFunctorName(const string& baseClassName1, const string& baseClassName2, const string& libName, shared_ptr<Functor> eu);
void storeFunctorName(const string& baseClassName1, const string& baseClassName2, const string& baseClassName3, const string& libName, shared_ptr<Functor> eu);
=== modified file 'core/Functor.hpp'
--- core/Functor.hpp 2009-12-01 14:56:39 +0000
+++ core/Functor.hpp 2009-12-09 16:43:25 +0000
@@ -14,11 +14,14 @@
#include<yade/lib-multimethods/FunctorWrapper.hpp>
class TimingDeltas;
+class Scene;
class Functor : public Serializable
{
public: virtual vector<std::string> getFunctorTypes(){throw;}
shared_ptr<TimingDeltas> timingDeltas;
+ //! updated before every dispatch loop by the dispatcher; DO NOT ABUSE access to scene, except for getting global variables like scene->dt.
+ Scene* scene;
// label to be able to retrieve an engine unit by its label
string label;
REGISTER_CLASS_AND_BASE(Functor,Serializable);
=== modified file 'pkg/common/Engine/Dispatcher/BoundDispatcher.cpp'
--- pkg/common/Engine/Dispatcher/BoundDispatcher.cpp 2009-12-04 23:46:45 +0000
+++ pkg/common/Engine/Dispatcher/BoundDispatcher.cpp 2009-12-09 16:43:25 +0000
@@ -17,9 +17,10 @@
CREATE_LOGGER(BoundDispatcher);
-void BoundDispatcher::action(Scene* ncb)
+void BoundDispatcher::action(Scene*)
{
- shared_ptr<BodyContainer>& bodies = ncb->bodies;
+ updateScenePtr();
+ shared_ptr<BodyContainer>& bodies = scene->bodies;
const long numBodies=(long)bodies->size();
bool haveBins=(bool)velocityBins;
if(sweepDist!=0 && haveBins){ LOG_FATAL("Only one of sweepDist or velocityBins can used!"); abort(); }
@@ -50,7 +51,7 @@
aabb->min=aabb->center-aabb->halfSize; aabb->max=aabb->center+aabb->halfSize;
}
}
- if(ncb->state && ncb->bound && ncb->shape) operator()(ncb->shape,ncb->bound,ncb->state->se3,ncb);
+ if(scene->state && scene->bound && scene->shape) operator()(scene->shape,scene->bound,scene->state->se3,scene);
}
=== modified file 'pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp'
--- pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2009-12-06 22:02:12 +0000
+++ pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2009-12-09 16:43:25 +0000
@@ -17,7 +17,7 @@
InteractionDispatchers::InteractionDispatchers(){
geomDispatcher=shared_ptr<InteractionGeometryDispatcher>(new InteractionGeometryDispatcher);
physDispatcher=shared_ptr<InteractionPhysicsDispatcher>(new InteractionPhysicsDispatcher);
- constLawDispatcher=shared_ptr<LawDispatcher>(new LawDispatcher);
+ lawDispatcher=shared_ptr<LawDispatcher>(new LawDispatcher);
alreadyWarnedNoCollider=false;
#ifdef IDISP_TIMING
timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
@@ -26,34 +26,37 @@
#define DISPATCH_CACHE
-void InteractionDispatchers::action(Scene* rootBody){
+void InteractionDispatchers::action(Scene*){
#ifdef IDISP_TIMING
timingDeltas->start();
#endif
- if(rootBody->interactions->unconditionalErasePending()>0 && !alreadyWarnedNoCollider){
+ if(scene->interactions->unconditionalErasePending()>0 && !alreadyWarnedNoCollider){
LOG_WARN("Interactions pending erase found (erased), no collider being used?");
alreadyWarnedNoCollider=true;
}
- Vector3r cellSize; if(rootBody->isPeriodic) cellSize=rootBody->cellMax-rootBody->cellMin;
- bool removeUnseenIntrs=(rootBody->interactions->iterColliderLastRun>=0 && rootBody->interactions->iterColliderLastRun==rootBody->currentIteration);
+ geomDispatcher->updateScenePtr();
+ physDispatcher->updateScenePtr();
+ lawDispatcher->updateScenePtr();
+ Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cellMax-scene->cellMin;
+ bool removeUnseenIntrs=(scene->interactions->iterColliderLastRun>=0 && scene->interactions->iterColliderLastRun==scene->currentIteration);
#ifdef YADE_OPENMP
- const long size=rootBody->interactions->size();
+ const long size=scene->interactions->size();
#pragma omp parallel for schedule(guided)
for(long i=0; i<size; i++){
- const shared_ptr<Interaction>& I=(*rootBody->interactions)[i];
+ const shared_ptr<Interaction>& I=(*scene->interactions)[i];
#else
- FOREACH(shared_ptr<Interaction> I, *rootBody->interactions){
+ FOREACH(shared_ptr<Interaction> I, *scene->interactions){
#endif
#ifdef DISPATCH_CACHE
- if(removeUnseenIntrs && !I->isReal() && I->iterLastSeen<rootBody->currentIteration) {
- rootBody->interactions->requestErase(I->getId1(),I->getId2());
+ if(removeUnseenIntrs && !I->isReal() && I->iterLastSeen<scene->currentIteration) {
+ scene->interactions->requestErase(I->getId1(),I->getId2());
continue;
}
- const shared_ptr<Body>& b1_=Body::byId(I->getId1(),rootBody);
- const shared_ptr<Body>& b2_=Body::byId(I->getId2(),rootBody);
+ const shared_ptr<Body>& b1_=Body::byId(I->getId1(),scene);
+ const shared_ptr<Body>& b2_=Body::byId(I->getId2(),scene);
- if(!b1_ || !b2_){ LOG_DEBUG("Body #"<<(b1_?I->getId2():I->getId1())<<" vanished, erasing intr #"<<I->getId1()<<"+#"<<I->getId2()<<"!"); rootBody->interactions->requestErase(I->getId1(),I->getId2(),/*force*/true); continue; }
+ if(!b1_ || !b2_){ LOG_DEBUG("Body #"<<(b1_?I->getId2():I->getId1())<<" vanished, erasing intr #"<<I->getId1()<<"+#"<<I->getId2()<<"!"); scene->interactions->requestErase(I->getId1(),I->getId2(),/*force*/true); continue; }
// go fast if this pair of bodies cannot interact at all
if((b1_->getGroupMask() & b2_->getGroupMask())==0) continue;
@@ -75,19 +78,19 @@
// and can call go in all cases
if(swap){I->swapOrder();}
// body pointers must be updated, in case we swapped
- const shared_ptr<Body>& b1=Body::byId(I->getId1(),rootBody);
- const shared_ptr<Body>& b2=Body::byId(I->getId2(),rootBody);
+ const shared_ptr<Body>& b1=Body::byId(I->getId1(),scene);
+ const shared_ptr<Body>& b2=Body::byId(I->getId2(),scene);
assert(I->functorCache.geom);
bool wasReal=I->isReal();
bool geomCreated;
- if(!rootBody->isPeriodic) geomCreated=I->functorCache.geom->go(b1->shape,b2->shape, *b1->state, *b2->state, Vector3r::ZERO, /*force*/false, I);
+ if(!scene->isPeriodic) geomCreated=I->functorCache.geom->go(b1->shape,b2->shape, *b1->state, *b2->state, Vector3r::ZERO, /*force*/false, I);
else{ // handle periodicity
Vector3r shift2(I->cellDist[0]*cellSize[0],I->cellDist[1]*cellSize[1],I->cellDist[2]*cellSize[2]);
geomCreated=I->functorCache.geom->go(b1->shape,b2->shape,*b1->state,*b2->state,shift2,/*force*/false,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
+ if(wasReal) scene->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
}
@@ -103,13 +106,13 @@
I->functorCache.phys->go(b1->material,b2->material,I);
assert(I->interactionPhysics);
- if(!wasReal) I->iterMadeReal=rootBody->currentIteration; // mark the interaction as created right now
+ if(!wasReal) I->iterMadeReal=scene->currentIteration; // mark the interaction as created right now
// LawDispatcher
// populating constLaw cache must be done after geom and physics dispatchers have been called, since otherwise the interaction
// would not have interactionGeometry and interactionPhysics yet.
if(!I->functorCache.constLaw){
- I->functorCache.constLaw=constLawDispatcher->getFunctor2D(I->interactionGeometry,I->interactionPhysics,swap);
+ I->functorCache.constLaw=lawDispatcher->getFunctor2D(I->interactionGeometry,I->interactionPhysics,swap);
if(!I->functorCache.constLaw){
LOG_FATAL("getFunctor2D returned empty functor for #"<<I->getId1()<<"+"<<I->getId2()<<", types "<<I->interactionGeometry->getClassName()<<"="<<I->interactionGeometry->getClassIndex()<<" and "<<I->interactionPhysics->getClassName()<<"="<<I->interactionPhysics->getClassIndex());
//abort();
@@ -118,26 +121,26 @@
assert(!swap); // reverse call would make no sense, as the arguments are of different types
}
assert(I->functorCache.constLaw);
- I->functorCache.constLaw->go(I->interactionGeometry,I->interactionPhysics,I.get(),rootBody);
+ I->functorCache.constLaw->go(I->interactionGeometry,I->interactionPhysics,I.get(),scene);
#else
- const shared_ptr<Body>& b1=Body::byId(I->getId1(),rootBody);
- const shared_ptr<Body>& b2=Body::byId(I->getId2(),rootBody);
+ const shared_ptr<Body>& b1=Body::byId(I->getId1(),scene);
+ const shared_ptr<Body>& b2=Body::byId(I->getId2(),scene);
// InteractionGeometryDispatcher
bool wasReal=I->isReal();
bool geomCreated =
b1->shape && b2->shape && // some bodies do not have shape
geomDispatcher->operator()(b1->shape, b2->shape, *b1->state, *b2->state, Vector3r::ZERO, I);
// FIXME: port from the part above
- if(rootBody->isPeriodic) { LOG_FATAL(__FILE__ ": Periodicity not handled without DISPATCH_CACHE."); abort(); }
+ if(scene->isPeriodic) { LOG_FATAL(__FILE__ ": Periodicity not handled without DISPATCH_CACHE."); abort(); }
if(!geomCreated){
- if(wasReal) *rootBody->interactions->requestErase(I->getId1(),I->getId2());
+ if(wasReal) *scene->interactions->requestErase(I->getId1(),I->getId2());
continue;
}
// InteractionPhysicsDispatcher
// geom may have swapped bodies, get bodies again
- physDispatcher->operator()(Body::byId(I->getId1(),rootBody)->material, Body::byId(I->getId2(),rootBody)->material,I);
+ physDispatcher->operator()(Body::byId(I->getId1(),scene)->material, Body::byId(I->getId2(),scene)->material,I);
// LawDispatcher
- constLawDispatcher->operator()(I->interactionGeometry,I->interactionPhysics,I.get(),rootBody);
+ lawDispatcher->operator()(I->interactionGeometry,I->interactionPhysics,I.get(),scene);
#endif
}
}
=== modified file 'pkg/common/Engine/Dispatcher/InteractionDispatchers.hpp'
--- pkg/common/Engine/Dispatcher/InteractionDispatchers.hpp 2009-12-04 23:07:34 +0000
+++ pkg/common/Engine/Dispatcher/InteractionDispatchers.hpp 2009-12-09 16:43:25 +0000
@@ -12,12 +12,12 @@
virtual void action(Scene*);
shared_ptr<InteractionGeometryDispatcher> geomDispatcher;
shared_ptr<InteractionPhysicsDispatcher> physDispatcher;
- shared_ptr<LawDispatcher> constLawDispatcher;
+ shared_ptr<LawDispatcher> lawDispatcher;
REGISTER_CLASS_AND_BASE(InteractionDispatchers,StandAloneEngine);
REGISTER_ATTRIBUTES(StandAloneEngine,
(geomDispatcher)
(physDispatcher)
- (constLawDispatcher)
+ (lawDispatcher)
);
DECLARE_LOGGER;
};
=== modified file 'pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp'
--- pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp 2009-12-06 22:02:12 +0000
+++ pkg/common/Engine/Dispatcher/InteractionGeometryDispatcher.cpp 2009-12-09 16:43:25 +0000
@@ -25,6 +25,7 @@
*/
shared_ptr<Interaction> InteractionGeometryDispatcher::explicitAction(const shared_ptr<Body>& b1, const shared_ptr<Body>& b2, bool force){
+ updateScenePtr();
if(force){
assert(b1->shape && b2->shape);
shared_ptr<Interaction> i(new Interaction(b1->getId(),b2->getId()));
@@ -38,44 +39,45 @@
}
}
-void InteractionGeometryDispatcher::action(Scene* ncb){
+void InteractionGeometryDispatcher::action(Scene*){
// Erase interaction that were requested for erase, but not processed by the collider, if any (and warn once about that, as it is suspicious)
- if(ncb->interactions->unconditionalErasePending()>0 && !alreadyWarnedNoCollider){
+ if(scene->interactions->unconditionalErasePending()>0 && !alreadyWarnedNoCollider){
LOG_WARN("Interactions pending erase found, no collider being used?");
alreadyWarnedNoCollider=true;
}
+ updateScenePtr();
- shared_ptr<BodyContainer>& bodies = ncb->bodies;
- Vector3r cellSize; if(ncb->isPeriodic) cellSize=ncb->cellMax-ncb->cellMin;
- bool removeUnseenIntrs=(ncb->interactions->iterColliderLastRun>=0 && ncb->interactions->iterColliderLastRun==ncb->currentIteration);
+ shared_ptr<BodyContainer>& bodies = scene->bodies;
+ Vector3r cellSize; if(scene->isPeriodic) cellSize=scene->cellMax-scene->cellMin;
+ bool removeUnseenIntrs=(scene->interactions->iterColliderLastRun>=0 && scene->interactions->iterColliderLastRun==scene->currentIteration);
#ifdef YADE_OPENMP
- const long size=ncb->interactions->size();
+ const long size=scene->interactions->size();
#pragma omp parallel for
for(long i=0; i<size; i++){
- const shared_ptr<Interaction>& I=(*ncb->interactions)[i];
+ const shared_ptr<Interaction>& I=(*scene->interactions)[i];
#else
- FOREACH(const shared_ptr<Interaction>& I, *ncb->interactions){
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
#endif
- if(removeUnseenIntrs && !I->isReal() && I->iterLastSeen<ncb->currentIteration) {
- ncb->interactions->requestErase(I->getId1(),I->getId2());
+ if(removeUnseenIntrs && !I->isReal() && I->iterLastSeen<scene->currentIteration) {
+ scene->interactions->requestErase(I->getId1(),I->getId2());
continue;
}
const shared_ptr<Body>& b1=(*bodies)[I->getId1()];
const shared_ptr<Body>& b2=(*bodies)[I->getId2()];
- if(!b1 || !b2){ LOG_DEBUG("Body #"<<(b1?I->getId2():I->getId1())<<" vanished, erasing intr #"<<I->getId1()<<"+#"<<I->getId2()<<"!"); ncb->interactions->requestErase(I->getId1(),I->getId2(),/*force*/true); continue; }
+ if(!b1 || !b2){ LOG_DEBUG("Body #"<<(b1?I->getId2():I->getId1())<<" vanished, erasing intr #"<<I->getId1()<<"+#"<<I->getId2()<<"!"); scene->interactions->requestErase(I->getId1(),I->getId2(),/*force*/true); continue; }
bool wasReal=I->isReal();
if (!b1->shape || !b2->shape) { assert(!wasReal); continue; } // some bodies do not have shape
bool geomCreated;
- if(!ncb->isPeriodic){
+ if(!scene->isPeriodic){
geomCreated=operator()(b1->shape, b2->shape, *b1->state, *b2->state, Vector3r::ZERO, /*force*/ false, I);
} else{
Vector3r shift2(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->shape, b2->shape, *b1->state, *b2->state, shift2, /*force*/ false, I);
}
// reset && erase interaction that existed but now has no geometry anymore
- if(wasReal && !geomCreated){ ncb->interactions->requestErase(I->getId1(),I->getId2()); }
+ if(wasReal && !geomCreated){ scene->interactions->requestErase(I->getId1(),I->getId2()); }
}
}
=== modified file 'pkg/common/Engine/Dispatcher/InteractionPhysicsDispatcher.cpp'
--- pkg/common/Engine/Dispatcher/InteractionPhysicsDispatcher.cpp 2009-12-06 22:02:12 +0000
+++ pkg/common/Engine/Dispatcher/InteractionPhysicsDispatcher.cpp 2009-12-09 16:43:25 +0000
@@ -17,21 +17,23 @@
void InteractionPhysicsDispatcher::explicitAction(shared_ptr<Material>& pp1, shared_ptr<Material>& pp2, shared_ptr<Interaction>& i){
// should we throw instead of asserting?
//assert(i->isReal());
+ updateScenePtr();
if(!i->interactionGeometry) throw runtime_error(string(__FILE__)+": explicitAction received interaction without interactionGeometry.");
operator()(pp1,pp2,i);
}
-void InteractionPhysicsDispatcher::action(Scene* ncb)
+void InteractionPhysicsDispatcher::action(Scene*)
{
- shared_ptr<BodyContainer>& bodies = ncb->bodies;
+ updateScenePtr();
+ shared_ptr<BodyContainer>& bodies = scene->bodies;
#ifdef YADE_OPENMP
- const long size=ncb->interactions->size();
+ const long size=scene->interactions->size();
#pragma omp parallel for
for(long i=0; i<size; i++){
- const shared_ptr<Interaction>& interaction=(*ncb->interactions)[i];
+ const shared_ptr<Interaction>& interaction=(*scene->interactions)[i];
#else
- FOREACH(const shared_ptr<Interaction>& interaction, *ncb->interactions){
+ FOREACH(const shared_ptr<Interaction>& interaction, *scene->interactions){
#endif
if(interaction->interactionGeometry){
shared_ptr<Body>& b1 = (*bodies)[interaction->getId1()];
@@ -39,7 +41,7 @@
bool hadPhys=interaction->interactionPhysics;
operator()(b1->material, b2->material, interaction);
assert(interaction->interactionPhysics);
- if(!hadPhys) interaction->iterMadeReal=ncb->currentIteration;
+ if(!hadPhys) interaction->iterMadeReal=scene->currentIteration;
}
}
}
=== modified file 'pkg/common/Engine/Dispatcher/LawDispatcher.cpp'
--- pkg/common/Engine/Dispatcher/LawDispatcher.cpp 2009-12-04 23:07:34 +0000
+++ pkg/common/Engine/Dispatcher/LawDispatcher.cpp 2009-12-09 16:43:25 +0000
@@ -1,18 +1,19 @@
// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
#include "LawDispatcher.hpp"
YADE_PLUGIN((LawDispatcher));
-void LawDispatcher::action(Scene* rootBody){
+void LawDispatcher::action(Scene*){
+ updateScenePtr();
#ifdef YADE_OPENMP
- const long size=rootBody->interactions->size();
+ const long size=scene->interactions->size();
#pragma omp parallel for
for(long i=0; i<size; i++){
- const shared_ptr<Interaction>& I=(*rootBody->interactions)[i];
+ const shared_ptr<Interaction>& I=(*scene->interactions)[i];
#else
- FOREACH(shared_ptr<Interaction> I, *rootBody->interactions){
+ FOREACH(shared_ptr<Interaction> I, *scene->interactions){
#endif
if(I->isReal()){
assert(I->interactionGeometry); assert(I->interactionPhysics);
- operator()(I->interactionGeometry,I->interactionPhysics,I.get(),rootBody);
+ operator()(I->interactionGeometry,I->interactionPhysics,I.get(),scene);
}
}
}
=== modified file 'pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp'
--- pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp 2009-12-05 23:12:46 +0000
+++ pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp 2009-12-09 16:43:25 +0000
@@ -59,18 +59,18 @@
#endif
}
-void NewtonsDampedLaw::action(Scene* scene)
+void NewtonsDampedLaw::action(Scene*)
{
scene->bex.sync();
Real dt=scene->dt;
// account for motion of the periodic boundary, if we remember its last position
// its velocity will count as max velocity of bodies
// otherwise the collider might not run if only the cell were changing without any particle motion
- if(scene->isPeriodic && prevCellSize!=Vector3r::ZERO){ maxVelocitySq=(prevCellSize-(scene->cellMax-scene->cellMin)).SquaredLength()/pow(dt,2); }
- else maxVelocitySq=0;
+ if(scene->isPeriodic && (prevCellMax!=scene->cellMax || prevCellMin!=scene->cellMin)){ cellChanged=true; maxVelocitySq=(prevCellMax-prevCellMin-(scene->cellMax-scene->cellMin)).SquaredLength()/pow(dt,2); }
+ else { maxVelocitySq=0; cellChanged=false; }
+ prevCellSize=prevCellMax-prevCellMin; cellSize=scene->cellMax-scene->cellMin;
haveBins=(bool)velocityBins;
if(haveBins) velocityBins->binVelSqInitialize(maxVelocitySq);
- if(scene->isPeriodic) prevCellSize=scene->cellMax-scene->cellMin;
#ifdef YADE_OPENMP
FOREACH(Real& thrMaxVSq, threadMaxVelocitySq) { thrMaxVSq=0; }
@@ -98,12 +98,12 @@
cundallDamp(dt,f,state->vel,state->accel);
leapfrogTranslate(scene,state,id,dt);
// rotate equation
- if (/*b->isSpheral || */ !exactAsphericalRot){ // spheral body or exactAsphericalRot disabled
+ if (!exactAsphericalRot){ // exactAsphericalRot disabled
const Vector3r& m=scene->bex.getTorque(id);
state->angAccel=diagDiv(m,state->inertia);
cundallDamp(dt,m,state->angVel,state->angAccel);
leapfrogSphericalRotate(scene,state,id,dt);
- } else { // non spheral body and exactAsphericalRot enabled
+ } else { // exactAsphericalRot enabled
const Vector3r& m=scene->bex.getTorque(id);
leapfrogAsphericalRotate(scene,state,id,dt,m);
}
@@ -152,13 +152,21 @@
FOREACH(const Real& thrMaxVSq, threadMaxVelocitySq) { maxVelocitySq=max(maxVelocitySq,thrMaxVSq); }
#endif
if(haveBins) velocityBins->binVelSqFinalize();
+ if(scene->isPeriodic) { prevCellMax=scene->cellMax; prevCellMin=scene->cellMin; }
}
inline void NewtonsDampedLaw::leapfrogTranslate(Scene* scene, State* state, const body_id_t& id, const Real& dt )
{
+ // for homothetic resize of the cell (if enabled), compute the position difference of the homothetic transformation
+ // the term ξ=(x'-xâ‚€')/(xâ‚'-xâ‚€') is normalized cell coordinate (' meaning at previous step),
+ // which is then used to compute new position x=ξ(xâ‚-xâ‚€)+xâ‚€. (per component)
+ // Then we update either velocity by (x-x')/Δt (homotheticCellResize==1)
+ // or position by (x-x') (homotheticCellResize==2)
+ Vector3r dPos;
+ if(cellChanged && homotheticCellResize){ for(int i=0; i<3; i++) dPos[i]=((state->pos[i]-prevCellMin[i])/(prevCellSize[i]))*cellSize[i]+scene->cellMin[i]-state->pos[i]; }
blockTranslateDOFs(state->blockedDOFs, state->accel);
- state->vel+=dt*state->accel;
- state->pos += state->vel*dt + scene->bex.getMove(id);
+ state->vel+=dt*state->accel; if(homotheticCellResize==1) state->vel+=dPos/dt;
+ state->pos += state->vel*dt + scene->bex.getMove(id); if(homotheticCellResize==2) state->pos+=dPos;
}
inline void NewtonsDampedLaw::leapfrogSphericalRotate(Scene* scene, State* state, const body_id_t& id, const Real& dt )
{
=== modified file 'pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp'
--- pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp 2009-12-05 23:12:46 +0000
+++ pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp 2009-12-09 16:43:25 +0000
@@ -50,7 +50,12 @@
Quaternionr DotQ(const Vector3r& angVel, const Quaternionr& Q);
inline void blockTranslateDOFs(unsigned blockedDOFs, Vector3r& v);
inline void blockRotateDOFs(unsigned blockedDOFs, Vector3r& v);
- Vector3r prevCellSize;
+ // cell corners from previous step, used to detect change, find max velocity and update positions if linearCellResize enabled
+ Vector3r prevCellMax, prevCellMin;
+ // cache sizes, recomputed before the loop at every step
+ Vector3r prevCellSize, cellSize;
+ // whether the cell has changed from the previous step
+ bool cellChanged;
@@ -61,6 +66,8 @@
Real maxVelocitySq;
/// Enable of the exact aspherical body rotation integrator
bool exactAsphericalRot;
+ //! Enable artificially moving all bodies with the periodic cell, such that its resizes are isotropic. 0: disabled (default), 1: position update, 2: velocity update.
+ int homotheticCellResize;
#ifdef YADE_OPENMP
vector<Real> threadMaxVelocitySq;
@@ -68,7 +75,7 @@
/// velocity bins (not used if not created)
shared_ptr<VelocityBins> velocityBins;
virtual void action(Scene *);
- NewtonsDampedLaw(): prevCellSize(Vector3r::ZERO),damping(0.2), maxVelocitySq(-1), exactAsphericalRot(false){
+ NewtonsDampedLaw(): prevCellSize(Vector3r::ZERO),damping(0.2), maxVelocitySq(-1), exactAsphericalRot(false), homotheticCellResize(0){
#ifdef YADE_OPENMP
threadMaxVelocitySq.resize(omp_get_max_threads());
#endif
=== modified file 'pkg/dem/PreProcessor/TriaxialTest.cpp'
--- pkg/dem/PreProcessor/TriaxialTest.cpp 2009-12-04 23:46:45 +0000
+++ pkg/dem/PreProcessor/TriaxialTest.cpp 2009-12-09 16:43:25 +0000
@@ -578,12 +578,12 @@
shared_ptr<InteractionDispatchers> ids(new InteractionDispatchers);
ids->geomDispatcher=interactionGeometryDispatcher;
ids->physDispatcher=interactionPhysicsDispatcher;
- ids->constLawDispatcher=shared_ptr<LawDispatcher>(new LawDispatcher);
+ ids->lawDispatcher=shared_ptr<LawDispatcher>(new LawDispatcher);
if(!facetWalls && !wallWalls){
shared_ptr<ef2_Spheres_Elastic_ElasticLaw> see(new ef2_Spheres_Elastic_ElasticLaw); see->sdecGroupMask=2;
- ids->constLawDispatcher->add(see);
+ ids->lawDispatcher->add(see);
} else {
- ids->constLawDispatcher->add(shared_ptr<Law2_Dem3Dof_Elastic_Elastic>(new Law2_Dem3Dof_Elastic_Elastic));
+ ids->lawDispatcher->add(shared_ptr<Law2_Dem3Dof_Elastic_Elastic>(new Law2_Dem3Dof_Elastic_Elastic));
}
rootBody->engines.push_back(ids);
} else {
=== modified file 'py/timing.py'
--- py/timing.py 2009-12-04 04:53:02 +0000
+++ py/timing.py 2009-12-09 16:43:25 +0000
@@ -60,7 +60,7 @@
if e.__class__.__name__=='InteractionDispatcher':
lines+=_engines_stats(e.geomDispatcher.functors,e.execTime,level+1)
lines+=_engines_stats(e.physDispatcher.functors,e.execTime,level+1)
- lines+=_engines_stats(e.constLawDispatcher.functors,e.execTime,level+1)
+ lines+=_engines_stats(e.lawDispatcher.functors,e.execTime,level+1)
elif e.__class__.__name__=='ParallelEngine': lines+=_engines_stats(e.slave,e.execTime,level+1)
if hereLines>1:
print _formatLine('TOTAL',totalTime,-1,totalTime,level); lines+=1
=== modified file 'py/yadeWrapper/yadeWrapper.cpp'
--- py/yadeWrapper/yadeWrapper.cpp 2009-12-08 12:08:48 +0000
+++ py/yadeWrapper/yadeWrapper.cpp 2009-12-09 16:43:25 +0000
@@ -344,7 +344,7 @@
list<shared_ptr<Functor> > eus;
FOREACH(const shared_ptr<Functor>& eu,ee->geomDispatcher->functorArguments) eus.push_back(eu);
FOREACH(const shared_ptr<Functor>& eu,ee->physDispatcher->functorArguments) eus.push_back(eu);
- FOREACH(const shared_ptr<Functor>& eu,ee->constLawDispatcher->functorArguments) eus.push_back(eu);
+ FOREACH(const shared_ptr<Functor>& eu,ee->lawDispatcher->functorArguments) eus.push_back(eu);
FOREACH(const shared_ptr<Functor>& eu,eus){
if(!eu->label.empty()){
PyGILState_STATE gstate; gstate = PyGILState_Ensure(); PyRun_SimpleString(("__builtins__."+eu->label+"=Omega().labeledEngine('"+eu->label+"')").c_str()); PyGILState_Release(gstate);
@@ -457,7 +457,7 @@
list<shared_ptr<Functor> > eus;
FOREACH(const shared_ptr<Functor>& eu,ee->geomDispatcher->functorArguments) eus.push_back(eu);
FOREACH(const shared_ptr<Functor>& eu,ee->physDispatcher->functorArguments) eus.push_back(eu);
- FOREACH(const shared_ptr<Functor>& eu,ee->constLawDispatcher->functorArguments) eus.push_back(eu);
+ FOREACH(const shared_ptr<Functor>& eu,ee->lawDispatcher->functorArguments) eus.push_back(eu);
FOREACH(const shared_ptr<Functor>& eu,eus){
if(eu->label==label) return python::object(eu);
}
@@ -622,7 +622,7 @@
shared_ptr<InteractionDispatchers> instance(new InteractionDispatchers);
FOREACH(shared_ptr<InteractionGeometryFunctor> gf, gff) instance->geomDispatcher->add(gf);
FOREACH(shared_ptr<InteractionPhysicsFunctor> pf, pff) instance->physDispatcher->add(pf);
- FOREACH(shared_ptr<LawFunctor> cf, cff) instance->constLawDispatcher->add(cf);
+ FOREACH(shared_ptr<LawFunctor> cf, cff) instance->lawDispatcher->add(cf);
return instance;
}
@@ -824,7 +824,7 @@
.def("__init__",python::make_constructor(InteractionDispatchers_ctor_lists))
.def_readonly("geomDispatcher",&InteractionDispatchers::geomDispatcher)
.def_readonly("physDispatcher",&InteractionDispatchers::physDispatcher)
- .def_readonly("constLawDispatcher",&InteractionDispatchers::constLawDispatcher);
+ .def_readonly("lawDispatcher",&InteractionDispatchers::lawDispatcher);
python::class_<ParallelEngine,shared_ptr<ParallelEngine>, python::bases<Engine>, noncopyable>("ParallelEngine")
.def("__init__",python::make_constructor(ParallelEngine_ctor_list))
.add_property("slaves",&ParallelEngine_slaves_get,&ParallelEngine_slaves_set);
=== modified file 'scripts/test-sphere-facet.py'
--- scripts/test-sphere-facet.py 2009-12-04 23:27:26 +0000
+++ scripts/test-sphere-facet.py 2009-12-09 16:43:25 +0000
@@ -42,13 +42,13 @@
O.dt=1e-4
print '** virgin dispatch matrix:'
-O.engines[3].constLawDispatcher.dump()
+O.engines[3].lawDispatcher.dump()
print '** class indices'
for c in 'Dem3DofGeom','Dem3DofGeom_FacetSphere','Dem3DofGeom_SphereSphere':
print eval(c)().classIndex,c
O.run(1000,True)
print '** used dispatch matrix'
-O.engines[3].constLawDispatcher.dump()
+O.engines[3].lawDispatcher.dump()
def setGravity():
