yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-03-01 14:29:51 +0100 (Sun, 01 Mar 2009)
New Revision: 1703

Modified:
   trunk/core/Interaction.cpp
   trunk/core/Interaction.hpp
   trunk/examples/collider-perf/perf.py
   trunk/extra/clump/Shop.cpp
   trunk/gui/py/utils.py
   trunk/gui/py/yade-multi
   trunk/lib/multimethods/DynLibDispatcher.hpp
   trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp
   trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp
   trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp
   trunk/pkg/dem/PreProcessor/TriaxialTest.cpp
Log:
1. Fix collider-perf (TriaxialTest has hardcoded dt=0.001 !!!!), set timestep by hand at first
2. Fix Shop:: loading spheres from file (regression from yesterday)
3. Add functorCache to Interaction and to InteractionDispatchers. About 5% improvement, but not yet finished.
4. Add DynLibDispatcher::getFunctor2D
5. Add BexContainer::sync() to other places in triaxial (thrown otherwise)




Modified: trunk/core/Interaction.cpp
===================================================================
--- trunk/core/Interaction.cpp	2009-03-01 01:41:37 UTC (rev 1702)
+++ trunk/core/Interaction.cpp	2009-03-01 13:29:51 UTC (rev 1703)
@@ -29,6 +29,7 @@
 	isReal = false;
 	isNeighbor = true;//NOTE : TriangulationCollider needs that
 
+	functorCache.geomExists=true;
 }
 
 

Modified: trunk/core/Interaction.hpp
===================================================================
--- trunk/core/Interaction.hpp	2009-03-01 01:41:37 UTC (rev 1702)
+++ trunk/core/Interaction.hpp	2009-03-01 13:29:51 UTC (rev 1703)
@@ -15,6 +15,8 @@
 typedef int body_id_t;
 
 class InteractionGeometryEngineUnit;
+class InteractionPhysicsEngineUnit;
+class ConstitutiveLaw;
 
 class Interaction : public Serializable
 {
@@ -39,6 +41,18 @@
 		//! swaps order of bodies within the interaction
 		void swapOrder();
 
+		//! cache functors that are called for this interaction. Currently used by InteractionDispatchers.
+		struct {
+			// Whether geometry dispatcher exists at all; this is different from !geom, since that can mean we haven't populated the cache yet.
+			// Therefore, geomExists must be initialized to true first (done in Interaction ctor).
+			bool geomExists;
+			// shared_ptr's are initialized to NULLs automagically
+			shared_ptr<InteractionGeometryEngineUnit> geom;
+			shared_ptr<InteractionPhysicsEngineUnit> phys;
+			shared_ptr<ConstitutiveLaw> constLaw;
+		} functorCache;
+			
+
 		#if 0
 			//! Whether both bodies involved in interaction satisfies given mask; provide rootBody for faster lookup
 			bool maskBothOK(int mask, MetaBody* rootBody=NULL){return (mask==0) || (Body::byId(id1,rootBody)->maskOK(mask) && Body::byId(id2,rootBody)->maskOK(mask));}

Modified: trunk/examples/collider-perf/perf.py
===================================================================
--- trunk/examples/collider-perf/perf.py	2009-03-01 01:41:37 UTC (rev 1702)
+++ trunk/examples/collider-perf/perf.py	2009-03-01 13:29:51 UTC (rev 1703)
@@ -17,6 +17,7 @@
 O.timingEnabled=True
 
 p=Preprocessor('TriaxialTest',{'importFilename':spheresFile}).load()
+O.dt=utils.PWaveTimeStep()
 utils.replaceCollider(StandAloneEngine(collider))
 
 O.step()

Modified: trunk/extra/clump/Shop.cpp
===================================================================
--- trunk/extra/clump/Shop.cpp	2009-03-01 01:41:37 UTC (rev 1702)
+++ trunk/extra/clump/Shop.cpp	2009-03-01 13:29:51 UTC (rev 1703)
@@ -485,7 +485,7 @@
 		int i=0;
 		FOREACH(const string& s, toks){
 			if(i<3) C[i]=lexical_cast<Real>(s);
-			if(i==4) r=lexical_cast<Real>(s);
+			if(i==3) r=lexical_cast<Real>(s);
 			i++;
 		}
 		if(i==0) continue; // empty line, skipped (can be the last one)

Modified: trunk/gui/py/utils.py
===================================================================
--- trunk/gui/py/utils.py	2009-03-01 01:41:37 UTC (rev 1702)
+++ trunk/gui/py/utils.py	2009-03-01 13:29:51 UTC (rev 1703)
@@ -354,7 +354,6 @@
 		for i in range(len(names)):
 			if names[i]=='description': o.tags['description']=values[i]
 			else:
-				print 'Parameter name:',names[i],names[i][0]
 				if names[i] not in kw.keys():
 					if (not unknownOk) and names[i][0]!='!': raise NameError("Parameter `%s' has no default value assigned"%names[i])
 				else: kw.pop(names[i])

Modified: trunk/gui/py/yade-multi
===================================================================
--- trunk/gui/py/yade-multi	2009-03-01 01:41:37 UTC (rev 1702)
+++ trunk/gui/py/yade-multi	2009-03-01 13:29:51 UTC (rev 1703)
@@ -14,7 +14,7 @@
 	def saveInfo(self):
 		log=file(self.log,'a')
 		log.write("""
-=================== JOB SUMMARY ================
+=================== JOB SUMMARY ================
 id      : %s
 status  : %d (%s)
 duration: %s

Modified: trunk/lib/multimethods/DynLibDispatcher.hpp
===================================================================
--- trunk/lib/multimethods/DynLibDispatcher.hpp	2009-03-01 01:41:37 UTC (rev 1702)
+++ trunk/lib/multimethods/DynLibDispatcher.hpp	2009-03-01 13:29:51 UTC (rev 1703)
@@ -442,7 +442,21 @@
 			#endif
 		  }
 		
-	public  : bool locateMultivirtualFunctor2D(int& index1, int& index2, shared_ptr<BaseClass1>& base1,shared_ptr<BaseClass2>& base2)
+	public:
+		/* Return pointer to the functor for two base classes given. Swap is true if the dispatch objects
+		 * should be swapped before calling Executor::go.
+		 */
+		shared_ptr<Executor> getFunctor2D(shared_ptr<BaseClass1>& base1, shared_ptr<BaseClass2>& base2, bool& swap){
+			int ix1, ix2;
+			if(!locateMultivirtualFunctor2D(ix1,ix2,base1,base2)){
+				return shared_ptr<Executor>();
+			}
+			swap=(bool)(callBacksInfo[ix1][ix2]);
+			return callBacks[ix1][ix2];
+		}
+	
+	
+		bool locateMultivirtualFunctor2D(int& index1, int& index2, shared_ptr<BaseClass1>& base1,shared_ptr<BaseClass2>& base2)
 		  {
 			index1 = base1->getClassIndex();
 			index2 = base2->getClassIndex();

Modified: trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp	2009-03-01 01:41:37 UTC (rev 1702)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp	2009-03-01 13:29:51 UTC (rev 1703)
@@ -8,6 +8,8 @@
 	constLawDispatcher=shared_ptr<ConstitutiveLawDispatcher>(new ConstitutiveLawDispatcher);
 }
 
+#define DISPATCH_CACHE
+
 void InteractionDispatchers::action(MetaBody* rootBody){
 	#ifdef YADE_OPENMP
 		const long size=rootBody->interactions->size();
@@ -17,16 +19,60 @@
 	#else
 		FOREACH(shared_ptr<Interaction> I, *rootBody->interactions){
 	#endif
+		#ifdef DISPATCH_CACHE
+			shared_ptr<Body>& b1=(*rootBody->bodies)[I->getId1()];
+			shared_ptr<Body>& b2=(*rootBody->bodies)[I->getId2()];
+			// we know there is no geometry functor already, take the short path
+			if(!I->functorCache.geomExists) { I->isReal=false; continue; }
+			// no interaction geometry for either of bodies; no interaction possible
+			if(!b1->interactingGeometry || !b2->interactingGeometry) { I->isReal=false; continue; }
+
 			// InteractionGeometryMetaEngine
+			if(!I->functorCache.geom || !I->functorCache.phys){
+				bool swap=false;
+				I->functorCache.geom=geomDispatcher->getFunctor2D(b1->interactingGeometry,b2->interactingGeometry,swap);
+				// returns NULL ptr if no functor exists; remember that and shortcut
+				if(!I->functorCache.geom) { I->functorCache.geomExists=false; continue; }
+				// arguments for the geom functor are in the reverse order (dispatcher would normally call goReverse).
+				// we don't remember the fact that is reverse, so we swap bodies within the interaction
+				// and can call go in all cases
+				if(swap) { I->swapOrder(); b1=(*rootBody->bodies)[I->getId1()]; b2=(*rootBody->bodies)[I->getId2()]; }
+			}
+			assert(I->functorCache.geom);
+			I->isReal=I->functorCache.geom->go(b1->interactingGeometry,b2->interactingGeometry,b1->physicalParameters->se3, b2->physicalParameters->se3,I);
+			if(!I->isReal) continue;
+
+			// InteractionPhysicsMetaEngine
+			if(!I->functorCache.phys){
+				bool swap=false; I->functorCache.phys=physDispatcher->getFunctor2D(b1->physicalParameters,b2->physicalParameters,swap);
+				assert(!swap); // InteractionPhysicsEngineUnits are symmetric
+			}
+			assert(I->functorCache.phys);
+			I->functorCache.phys->go(b1->physicalParameters,b2->physicalParameters,I);
+
+			// ConstitutiveLawDispatcher
+			// 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){
+				bool swap=false; I->functorCache.constLaw=constLawDispatcher->getFunctor2D(I->interactionGeometry,I->interactionPhysics,swap);
+				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);
+
+		#else
 			const shared_ptr<Body>& b1=Body::byId(I->getId1(),rootBody);
 			const shared_ptr<Body>& b2=Body::byId(I->getId2(),rootBody);
+			// InteractionGeometryMetaEngine
 			I->isReal =
 				b1->interactingGeometry && b2->interactingGeometry && // some bodies do not have interactingGeometry
 				geomDispatcher->operator()(b1->interactingGeometry, b2->interactingGeometry, b1->physicalParameters->se3, b2->physicalParameters->se3,I);
 			if(!I->isReal) continue;
 			// InteractionPhysicsMetaEngine
-			physDispatcher->operator()(b1->physicalParameters, b2->physicalParameters,I);
+			// geom may have swapped bodies, get bodies again
+			physDispatcher->operator()(Body::byId(I->getId1(),rootBody)->physicalParameters, Body::byId(I->getId2(),rootBody)->physicalParameters,I);
 			// ConstitutiveLawDispatcher
 			constLawDispatcher->operator()(I->interactionGeometry,I->interactionPhysics,I.get(),rootBody);
+		#endif
 		}
 }

Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp	2009-03-01 01:41:37 UTC (rev 1702)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp	2009-03-01 13:29:51 UTC (rev 1703)
@@ -400,6 +400,9 @@
  */
 Real TriaxialStressController::ComputeUnbalancedForce(MetaBody * ncb, bool maxUnbalanced)
 {
+	#ifdef BEX_CONTAINER
+		ncb->bex.sync();
+	#endif
 	//compute the mean contact force
 	Real MeanForce = 0.f;
 	long nForce = 0;

Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp	2009-03-01 01:41:37 UTC (rev 1702)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp	2009-03-01 13:29:51 UTC (rev 1703)
@@ -104,7 +104,7 @@
 		void updateStiffness(MetaBody* ncb);
 		void computeStressStrain(MetaBody* ncb); //Compute stresses on walls as "Vector3r stress[6]", compute meanStress, strain[3] and mean strain
 		//! Compute the mean/max unbalanced force in the assembly (normalized by mean contact force)
-    		Real ComputeUnbalancedForce(MetaBody * ncb, bool maxUnbalanced=false);
+    	Real ComputeUnbalancedForce(MetaBody * ncb, bool maxUnbalanced=false);
 
 		DECLARE_LOGGER;
 		

Modified: trunk/pkg/dem/PreProcessor/TriaxialTest.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/TriaxialTest.cpp	2009-03-01 01:41:37 UTC (rev 1702)
+++ trunk/pkg/dem/PreProcessor/TriaxialTest.cpp	2009-03-01 13:29:51 UTC (rev 1703)
@@ -266,26 +266,32 @@
 	/* if _mean_radius is not given (i.e. <=0), then calculate it from box size;
 	 * OTOH, if it is specified, scale the box preserving its ratio and lowerCorner so that the radius can be as requested
 	 */
-	Vector3r dimensions=upperCorner-lowerCorner; Real volume=dimensions.X()*dimensions.Y()*dimensions.Z();
 	Real porosity=.75;
-	Real really_radiusMean;
+	
+	vector<BasicSphere> sphere_list;
 
-	if(radiusMean<=0) really_radiusMean=pow(volume*(1-porosity)/(Mathr::PI*(4/3.)*numberOfGrains),1/3.);
+	if(importFilename==""){
+		Vector3r dimensions=upperCorner-lowerCorner; Real volume=dimensions.X()*dimensions.Y()*dimensions.Z();
+		Real really_radiusMean;
+		if(radiusMean<=0) really_radiusMean=pow(volume*(1-porosity)/(Mathr::PI*(4/3.)*numberOfGrains),1/3.);
+		else {
+			bool fixedDims[3];
+			fixedDims[0]=fixedBoxDims.find('x')!=string::npos; fixedDims[1]=fixedBoxDims.find('y')!=string::npos; fixedDims[2]=fixedBoxDims.find('z')!=string::npos;
+			int nScaled=(3-(int)fixedDims[0]+(int)fixedDims[1]+(int)fixedDims[2]);
+			if(nScaled==0) throw std::invalid_argument("At most 2 (not 3) axes can have fixed dimensions in fixedBoxDims if scaling for given radiusMean.");
+			Real boxScaleFactor=radiusMean*pow((4/3.)*Mathr::PI*numberOfGrains/(volume*(1-porosity)),1./nScaled);
+			LOG_INFO("Mean radius value of "<<radiusMean<<" requested, scaling "<<nScaled<<" dimensions by "<<boxScaleFactor);
+			dimensions[0]*=fixedDims[0]?1.:boxScaleFactor; dimensions[1]*=fixedDims[1]?1.:boxScaleFactor; dimensions[2]*=fixedDims[2]?1.:boxScaleFactor;
+			upperCorner=lowerCorner+dimensions;
+			really_radiusMean=radiusMean;
+		}
+		message+=GenerateCloud(sphere_list, lowerCorner, upperCorner, numberOfGrains, radiusStdDev, really_radiusMean, porosity);
+	}
 	else {
-		bool fixedDims[3];
-		fixedDims[0]=fixedBoxDims.find('x')!=string::npos; fixedDims[1]=fixedBoxDims.find('y')!=string::npos; fixedDims[2]=fixedBoxDims.find('z')!=string::npos;
-		int nScaled=(3-(int)fixedDims[0]+(int)fixedDims[1]+(int)fixedDims[2]);
-		if(nScaled==0) throw std::invalid_argument("At most 2 (not 3) axes can have fixed dimensions in fixedBoxDims if scaling for given radiusMean.");
-		Real boxScaleFactor=radiusMean*pow((4/3.)*Mathr::PI*numberOfGrains/(volume*(1-porosity)),1./nScaled);
-		LOG_INFO("Mean radius value of "<<radiusMean<<" requested, scaling "<<nScaled<<" dimensions by "<<boxScaleFactor);
-		dimensions[0]*=fixedDims[0]?1.:boxScaleFactor; dimensions[1]*=fixedDims[1]?1.:boxScaleFactor; dimensions[2]*=fixedDims[2]?1.:boxScaleFactor;
-		upperCorner=lowerCorner+dimensions;
-		really_radiusMean=radiusMean;
+		if(radiusMean>0) LOG_WARN("radiusMean ignored, since importFilename specified.");
+		sphere_list=Shop::loadSpheresFromFile(importFilename,lowerCorner,upperCorner);
 	}
-
 	
-	
-	
 	if(boxWalls)
 	{
 	// bottom box
@@ -388,11 +394,7 @@
 			}
 			 
 	}
-	
-	
-	vector<BasicSphere> sphere_list;
-	if(importFilename!="") sphere_list=Shop::loadSpheresFromFile(importFilename,lowerCorner,upperCorner);
-	else message+=GenerateCloud(sphere_list, lowerCorner, upperCorner, numberOfGrains, radiusStdDev, really_radiusMean, porosity);
+
 	vector<BasicSphere>::iterator it = sphere_list.begin();
 	vector<BasicSphere>::iterator it_end = sphere_list.end();
 	FOREACH(const BasicSphere& it, sphere_list){
@@ -401,10 +403,8 @@
 		rootBody->bodies->insert(body);
 	}	
 
-	
-	
+
 	return true;
-
 }