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;
-
}
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
