1. Invalidate persistent collider data when interactions are cleared (not tested yet)

2. Add renderer pointer to scene, so that functors can find display attributes
3. Gl1_NormPhys honors displacement scaling now
4. Fix scalarOnColorScale bug
5. Add RadialForceEngine, fix AxialGravityEngine

core/InteractionContainer.cpp

@@ -33,6 +33,7 @@ void InteractionContainer::clear(){
 	linIntrs.clear(); // clear the linear container
 	pendingErase.clear();
 	currSize=0;
+	dirty=true;
 }
 
 

core/InteractionContainer.hpp

@@ -60,8 +60,10 @@ class InteractionContainer: public Serializable{
 		// used only during serialization/deserialization
 		vector<shared_ptr<Interaction> > interaction;
 	public:
+		// flag for notifying the collider that persistent data should be invalidated
+		bool dirty;
 		// required by the class factory... :-|
-		InteractionContainer(): currSize(0),serializeSorted(false),iterColliderLastRun(-1){
+		InteractionContainer(): currSize(0),dirty(false),serializeSorted(false),iterColliderLastRun(-1){
 			bodies=NULL;
 			#ifdef YADE_OPENMP
 				threadsPendingErase.resize(omp_get_max_threads());
@@ -159,7 +161,7 @@ class InteractionContainer: public Serializable{
 	void postSave(InteractionContainer&);
 
 
-	REGISTER_ATTRIBUTES(Serializable,(interaction)(serializeSorted));
+	REGISTER_ATTRIBUTES(Serializable,(interaction)(serializeSorted)(dirty));
 	REGISTER_CLASS_AND_BASE(InteractionContainer,Serializable);
 };
 REGISTER_SERIALIZABLE(InteractionContainer);

core/Scene.cpp

@@ -74,10 +74,14 @@ void Scene::moveToNextTimeStep(){
 		checkStateTypes();
 		forces.resize(bodies->size()); // optimization, not necessary
 	}
-	// substepping or not, update engines from _nextEngines, if defined
-	if(!_nextEngines.empty() && subStep<0){
+	// substepping or not, update engines from _nextEngines, if defined, at the beginning of step
+	// subStep can be 0, which happens if simulations is saved in the middle of step (without substepping)
+	// this assumes that prologue will not set _nextEngines, which is safe hopefully
+	if(!_nextEngines.empty() && (subStep<0 || (subStep<=0 && !subStepping))){
 		engines=_nextEngines;
 		_nextEngines.clear();
+		// hopefully this will not break in some margin cases (subStepping with setting _nextEngines and such)
+		subStep=-1;
 	}
 	if(likely(!subStepping && subStep<0)){
 		/* set substep to 0 during the loop, so that engines/nextEngines handler know whether we are inside the loop currently */

core/Scene.hpp

@@ -26,6 +26,9 @@
 
 
 class Bound;
+#ifdef YADE_OPENGL
+	class OpenGLRenderer;
+#endif
 
 class Scene: public Serializable{
 	public:
@@ -54,6 +57,10 @@ class Scene: public Serializable{
 
 		shared_ptr<Engine> engineByName(const string& s);
 
+		#ifdef YADE_OPENGL
+			shared_ptr<OpenGLRenderer> renderer;
+		#endif
+
 		void postLoad(Scene&);
 
 		// bits for Scene::flags

gui/qt4/GLViewer.cpp

@@ -515,13 +515,16 @@ void GLViewer::draw()
 			}
 			renderer->clipPlaneSe3[manipulatedClipPlane]=newSe3;
 		}
-		renderer->render(Omega::instance().getScene(), selectedName());
+		const shared_ptr<Scene>& scene=Omega::instance().getScene();
+		scene->renderer=renderer;
+		renderer->render(scene, selectedName());
 	}
 }
 
 void GLViewer::drawWithNames(){
 	qglviewer::Vec vd=camera()->viewDirection(); renderer->viewDirection=Vector3r(vd[0],vd[1],vd[2]);
 	const shared_ptr<Scene> scene(Omega::instance().getScene());
+	scene->renderer=renderer;
 	renderer->scene=scene;
 	renderer->renderShape();
 }

lib/base/Logging.hpp

@@ -61,8 +61,8 @@
 #	define _LOG_HEAD __FILE__ ":"<<__LINE__<<" "<<__FUNCTION__<<": "
 #	define LOG_TRACE(msg) // _POOR_MANS_LOG("TRACE",msg)
 #	define LOG_DEBUG(msg) // _POOR_MANS_LOG("DEBUG",msg)
-#	define LOG_INFO(msg) // _POOR_MANS_LOG("INFO ",msg)
-#	define LOG_WARN(msg) // _POOR_MANS_LOG("WARN ",msg)
+#	define LOG_INFO(msg)  // _POOR_MANS_LOG("INFO ",msg)
+#	define LOG_WARN(msg)  _POOR_MANS_LOG("WARN ",msg)
 #	define LOG_ERROR(msg) _POOR_MANS_LOG("ERROR",msg)
 #	define LOG_FATAL(msg) _POOR_MANS_LOG("FATAL",msg)
 

lib/serialization/Serializable.hpp

@@ -50,7 +50,7 @@ using namespace std;
 namespace yade{
 	namespace Attr{
 		// keep in sync with py/wrapper/yadeWrapper.cpp !
-		enum flags { noSave=1, readonly=2, triggerPostLoad=4, hidden=8, noResize=16, };
+		enum flags { noSave=1, readonly=2, triggerPostLoad=4, hidden=8, noResize=16 };
 	};
 };
 using namespace yade;

pkg/common/ForceEngine.cpp

@@ -10,7 +10,7 @@
 #include<yade/core/IGeom.hpp>
 #include<yade/core/IPhys.hpp>
 
-YADE_PLUGIN((ForceEngine)(InterpolatingDirectedForceEngine));
+YADE_PLUGIN((ForceEngine)(InterpolatingDirectedForceEngine)(RadialForceEngine));
 
 void ForceEngine::action(){
 	FOREACH(Body::id_t id, ids){
@@ -26,4 +26,14 @@ void InterpolatingDirectedForceEngine::action(){
 	ForceEngine::action();
 }
 
+void RadialForceEngine::postLoad(RadialForceEngine&){ axisDir.normalize(); }
 
+void RadialForceEngine::action(){
+	FOREACH(Body::id_t id, ids){
+		assert(scene->bodies->exists(id));
+		const Vector3r& pos=Body::byId(id,scene)->state->pos;
+		Vector3r radial=(pos - (axisPt+axisDir * /* t */ ((pos-axisPt).dot(axisDir)))).normalized();
+		if(radial.squaredNorm()==0) continue;
+		scene->forces.addForce(id,fNorm*radial);
+	}
+}

pkg/common/ForceEngine.hpp

@@ -38,4 +38,13 @@ class InterpolatingDirectedForceEngine: public ForceEngine{
 };
 REGISTER_SERIALIZABLE(InterpolatingDirectedForceEngine);
 
-
+struct RadialForceEngine: public PartialEngine{
+	virtual void action();
+	virtual void postLoad(RadialForceEngine&);
+	YADE_CLASS_BASE_DOC_ATTRS(RadialForceEngine,PartialEngine,"Apply force of given magnitude directed away from spatial axis.",
+		((Vector3r,axisPt,Vector3r::Zero(),,"Point on axis"))
+		((Vector3r,axisDir,Vector3r::UnitX(),Attr::triggerPostLoad,"Axis direction (normalized automatically)"))
+		((Real,fNorm,0,,"Applied force magnitude"))
+	);
+};
+REGISTER_SERIALIZABLE(RadialForceEngine);

pkg/common/GLDrawFunctors.hpp

@@ -22,6 +22,8 @@ struct GLViewInfo{
 	Real sceneRadius;
 };
 
+class OpenGLRenderer;
+
 #define GL_FUNCTOR(Klass,typelist,renderedType) class Klass: public Functor1D<renderedType,void,typelist>{public:\
 	virtual ~Klass(){};\
 	virtual string renders() const { throw std::runtime_error(#Klass ": unregistered gldraw class.\n"); };\

pkg/common/Gl1_NormPhys.cpp

@@ -2,10 +2,12 @@
 
 #include<yade/core/Scene.hpp>
 #include<yade/pkg/common/Gl1_NormPhys.hpp>
+#include<yade/pkg/common/OpenGLRenderer.hpp>
 #include<yade/pkg/common/NormShearPhys.hpp>
 #include<yade/pkg/dem/DemXDofGeom.hpp>
 #include<yade/pkg/dem/Shop.hpp>
 
+
 YADE_PLUGIN((Gl1_NormPhys));
 
 GLUquadric* Gl1_NormPhys::gluQuadric=NULL;
@@ -68,8 +70,15 @@ void Gl1_NormPhys::go(const shared_ptr<IPhys>& ip, const shared_ptr<Interaction>
 		// max(r,0) handles r<0 which is the case for "radius" of the facet in Dem3DofGeom_FacetSphere
 		Vector3r cp=scene->isPeriodic? scene->cell->wrapShearedPt(geom->contactPoint) : geom->contactPoint;
 		Vector3r p1=cp-max(geom->refR1,0.)*geom->normal;
-		Vector3r relPos=/*p2*/(cp+max(geom->refR2,0.)*geom->normal)-p1;
-		Real dist=max(geom->refR1,0.)+max(geom->refR2,0.);
+		Vector3r p2=cp+max(geom->refR2,0.)*geom->normal;
+		const Vector3r& dispScale=scene->renderer ? scene->renderer->dispScale : Vector3r::Ones(); 
+		if(dispScale!=Vector3r::Ones()){
+			// move p1 and p2 by the same amounts as particles themselves would be moved
+			p1+=dispScale.cwise()*Vector3r(b1->state->pos-b1->state->refPos);
+			p2+=dispScale.cwise()*Vector3r(b2->state->pos-b2->state->refPos);
+		}
+		Vector3r relPos=p2-p1;
+		Real dist=relPos.norm(); //max(geom->refR1,0.)+max(geom->refR2,0.);
 	#endif
 
 

pkg/common/GravityEngines.cpp

@@ -44,8 +44,9 @@ void AxialGravityEngine::action(){
 		const Vector3r& x0=b->state->pos;
 		const Vector3r& x1=axisPoint;
 		const Vector3r x2=axisPoint+axisDirection;
-		Vector3r closestAxisPoint=(x2-x1) * /* t */ (-(x1-x0).dot(x2-x1))/((x2-x1).squaredNorm());
+		Vector3r closestAxisPoint=x1+(x2-x1) * /* t */ (-(x1-x0).dot(x2-x1))/((x2-x1).squaredNorm());
 		Vector3r toAxis=closestAxisPoint-x0; toAxis.normalize();
+		if(toAxis.squaredNorm()==0) continue;
 		scene->forces.addForce(b->getId(),acceleration*b->state->mass*toAxis);
 	}
 }

pkg/common/InsertionSortCollider.cpp

@@ -106,6 +106,7 @@ vector<Body::id_t> InsertionSortCollider::probeBoundingVolume(const Bound& bv){
 			if(fastestBodyMaxDist>=verletDist) return true;
 		}
 		if((size_t)BB[0].size!=2*scene->bodies->size()) return true;
+		if(scene->interactions->dirty) 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
 		scene->interactions->erasePending(*this,scene);
@@ -162,6 +163,12 @@ void InsertionSortCollider::action(){
 		boundDispatcher->scene=scene;
 		boundDispatcher->action();
 
+		// if interactions are dirty, force reinitialization
+		if(scene->interactions->dirty){
+			doInitSort=true;
+			scene->interactions->dirty=false;
+		}
+
 		if(verletDist<0){
 			Real minR=std::numeric_limits<Real>::infinity();
 			FOREACH(const shared_ptr<Body>& b, *scene->bodies){

pkg/common/InteractionLoop.cpp

@@ -35,6 +35,9 @@ void InteractionLoop::action(){
 		LOG_WARN("Interactions pending erase found (erased), no collider being used?");
 		alreadyWarnedNoCollider=true;
 	}
+	if(scene->interactions->dirty){
+		throw std::logic_error("InteractionContainer::dirty is true; the collider should re-initialize in such case and clear the dirty flag.");
+	}
 	// update Scene* of the dispatchers
 	geomDispatcher->scene=physDispatcher->scene=lawDispatcher->scene=scene;
 	// ask dispatchers to update Scene* of their functors

pkg/dem/Shop.cpp

@@ -428,9 +428,9 @@ void Shop::getViscoelasticFromSpheresInteraction( Real tc, Real en, Real es, sha
  */
 Vector3r Shop::scalarOnColorScale(Real x, Real xmin, Real xmax){
 	Real xnorm=min((Real)1.,max((x-xmin)/(xmax-xmin),(Real)0.));
-	if(xnorm<.25) return Vector3r(0,.4*xnorm,1);
+	if(xnorm<.25) return Vector3r(0,4.*xnorm,1);
 	if(xnorm<.5)  return Vector3r(0,1,1.-4.*(xnorm-.25));
-	if(xnorm<.75) return Vector3r(4*(xnorm-.5),1.,0);
+	if(xnorm<.75) return Vector3r(4.*(xnorm-.5),1.,0);
 	return Vector3r(1,1-4*(xnorm-.75),0);
 }
 

py/wrapper/yadeWrapper.cpp

@@ -605,7 +605,7 @@ BOOST_PYTHON_MODULE(wrapper)
 		.def("eraseNonReal",&pyInteractionContainer::eraseNonReal,"Erase all interactions that are not :yref:`real <InteractionContainer.isReal>`.")
 		.def("erase",&pyInteractionContainer::erase,"Erase one interaction, given by id1, id2 (internally, ``requestErase`` is called -- the interaction might still exist as potential, if the :yref:`Collider` decides so).")
 		.add_property("serializeSorted",&pyInteractionContainer::serializeSorted_get,&pyInteractionContainer::serializeSorted_set)
-		.def("clear",&pyInteractionContainer::clear,"Remove all interactions");
+		.def("clear",&pyInteractionContainer::clear,"Remove all interactions, and invalidate persistent collider data (if the collider supports it).");
 	python::class_<pyInteractionIterator>("InteractionIterator",python::init<pyInteractionIterator&>())
 		.def("__iter__",&pyInteractionIterator::pyIter)
 		.def("next",&pyInteractionIterator::pyNext);

scripts/test/beam-l6geom.py

@@ -33,9 +33,9 @@ try:
 	rr=qt.Renderer()
 	rr.intrGeom=True
 	Gl1_L6Geom.phiScale=30; Gl1_L3Geom.uScale=20
-	O.engines=O.engines+[
-		qt.SnapshotEngine(fileBase=O.tmpFilename(),label='snapper',iterPeriod=300,deadTimeout=20),
-		PyRunner(iterPeriod=330000,command='utils.makeVideo(snapper.snapshots,out="beam-l6geom.avi"); snapper.dead=True; O.pause()')
-	]
+	#O.engines=O.engines+[
+	#	qt.SnapshotEngine(fileBase=O.tmpFilename(),label='snapper',iterPeriod=300,deadTimeout=20),
+	#	PyRunner(iterPeriod=330000,command='utils.makeVideo(snapper.snapshots,out="beam-l6geom.avi"); snapper.dead=True; O.pause()')
+	#]
 except ImportError: pass
 O.run()