yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01008
[svn] r1680 - in trunk: . core extra extra/clump extra/usct gui/py pkg/common/Engine/DeusExMachina pkg/common/Engine/MetaEngine pkg/common/Engine/StandAloneEngine pkg/dem/Engine/DeusExMachina pkg/dem/Engine/EngineUnit pkg/dem/Engine/StandAloneEngine
Author: eudoxos
Date: 2009-02-22 21:08:45 +0100 (Sun, 22 Feb 2009)
New Revision: 1680
Modified:
trunk/SConstruct
trunk/core/BexContainer.hpp
trunk/core/Omega.hpp
trunk/extra/Brefcom.cpp
trunk/extra/clump/Shop.cpp
trunk/extra/clump/Shop.hpp
trunk/extra/usct/UniaxialStrainControlledTest.cpp
trunk/gui/py/_utils.cpp
trunk/gui/py/yadeControl.cpp
trunk/pkg/common/Engine/DeusExMachina/GravityEngines.cpp
trunk/pkg/common/Engine/MetaEngine/ConstitutiveLaw.hpp
trunk/pkg/common/Engine/StandAloneEngine/PhysicalActionContainerInitializer.cpp
trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp
trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp
trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
Log:
1. Add BexContainer (same interface, but separate implementation) in openMP flavor.
2. openmp=0 is the default for scons now (for some time), since openMP BexContainer gives about 7% slowdown and thus is not useful before at least some engines are parallelized as well.
3. Removed getting references to Force/Torque from BexContainer, since that breaks concurrency. Instead, addForce/addTorque must be used. If you need to read back with getForce/getTorque, you have to call sync() (expensive!!) beforehand, which will compute summary values for all bodies; if you forget, exception will be thrown. Note that sync() is invalidated at next write operation.
4. Adapted a few constitutive laws and other engines for the changes mentiones here.
Modified: trunk/SConstruct
===================================================================
--- trunk/SConstruct 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/SConstruct 2009-02-22 20:08:45 UTC (rev 1680)
@@ -102,8 +102,8 @@
print '@@@ Using profile',profile,'('+optsFile+') @@@'
# defaults for various profiles
-if profile=='default': defOptions={'debug':1,'variant':'','optimize':0,'openmp':True}
-elif profile=='opt': defOptions={'debug':0,'variant':'-opt','optimize':1,'openmp':True}
+if profile=='default': defOptions={'debug':1,'variant':'','optimize':0,'openmp':False}
+elif profile=='opt': defOptions={'debug':0,'variant':'-opt','optimize':1,'openmp':False}
else: defOptions={'debug':0,'optimize':0,'variant':profile,'openmp':True}
@@ -390,7 +390,7 @@
### COMPILER
if env['debug']: env.Append(CXXFLAGS='-ggdb3',CPPDEFINES=['YADE_DEBUG'])
else: env.Append(CXXFLAGS='-O2')
-if env['openmp']: env.Append(CXXFLAGS='-fopenmp',LIBS='gomp')
+if env['openmp']: env.Append(CXXFLAGS='-fopenmp',LIBS='gomp',CPPDEFINES='YADE_OPENMP')
if env['optimize']:
env.Append(CXXFLAGS=Split('-O3 -ffast-math -march=%s'%env['march']),
CPPDEFINES=[('YADE_CAST','static_cast'),('YADE_PTR_CAST','static_pointer_cast'),'NDEBUG'])
Modified: trunk/core/BexContainer.hpp
===================================================================
--- trunk/core/BexContainer.hpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/core/BexContainer.hpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -7,6 +7,83 @@
// for body_id_t
#include<yade/core/Interaction.hpp>
+#ifdef YADE_OPENMP
+
+#include<omp.h>
+
+class BexContainer{
+ private:
+ typedef std::vector<Vector3r> vvector;
+ std::vector<vvector> _forceData;
+ std::vector<vvector> _torqueData;
+ vvector _force, _torque;
+ size_t size;
+ int nThreads;
+ bool synced;
+ boost::mutex globalMutex;
+
+ inline void ensureSize(body_id_t id){
+ assert(nThreads>omp_get_thread_num());
+ if (size<=(size_t)id) resize(min((size_t)1.5*(id+100),(size_t)(id+2000)));
+ }
+
+ inline void ensureSynced(){ if(!synced) throw runtime_error("BexContainer not thread-synchronized; call sync() first!"); }
+ public:
+ BexContainer(): size(0), synced(true){
+ nThreads=omp_get_max_threads();
+ for(int i=0; i<nThreads; i++){
+ _forceData.push_back(vvector()); _torqueData.push_back(vvector());
+ }
+ }
+
+ /* To be benchmarked: sum thread data in getForce/getTorque upon request for each body individually instead of by the sync() function globally */
+ const Vector3r& getForce(body_id_t id) { ensureSize(id); ensureSynced(); return _force[id]; }
+ void addForce(body_id_t id, const Vector3r& f){ ensureSize(id); synced=false; _forceData[omp_get_thread_num()][id]+=f;}
+ const Vector3r& getTorque(body_id_t id) { ensureSize(id); ensureSynced(); return _torque[id]; }
+ void addTorque(body_id_t id, const Vector3r& f){ ensureSize(id); synced=false; _torqueData[omp_get_thread_num()][id]+=f;}
+
+ /* Sum contributions from all threads, save to the 0th thread storage.
+ * Locks globalMutex, since one thread modifies other threads' data.
+ * Must be called before get* methods are used. Exception is thrown otherwise, since data are not consistent.
+ */
+ inline void sync(){
+ if(synced) return;
+ boost::mutex::scoped_lock lock(globalMutex);
+ // #pragma omp parallel for private(sumF,sumT,thread);
+ for(size_t id=0; id<size; id++){
+ Vector3r sumF(Vector3r::ZERO), sumT(Vector3r::ZERO);
+ for(int thread=0; thread<omp_get_max_threads(); thread++){ sumF+=_forceData[thread][id]; sumT+=_torqueData[thread][id];}
+ _force[id]=sumF; _torque[id]=sumT;
+ }
+ synced=true;
+ }
+
+ /* Change size of containers (number of bodies).
+ * Locks globalMutex, since on threads modifies other threads' data.
+ * Called very rarely (a few times at the beginning of the simulation)
+ */
+ void resize(size_t newSize){
+ boost::mutex::scoped_lock lock(globalMutex);
+ if(size>=newSize) return; // in case on thread was waiting for resize, but it was already satisfied by another one
+ for(int thread=0; thread<nThreads; thread++){
+ _forceData [thread].resize(newSize);
+ _torqueData[thread].resize(newSize);
+ }
+ _force.resize(newSize); _torque.resize(newSize);
+ size=newSize;
+ }
+
+ void reset(){
+ for(int thread=0; thread<nThreads; thread++){
+ memset(_forceData [thread][0], 0,sizeof(Vector3r)*size);
+ memset(_torqueData[thread][0],0,sizeof(Vector3r)*size);
+ }
+ memset(_force [0], 0,sizeof(Vector3r)*size); memset(_torque[0], 0,sizeof(Vector3r)*size);
+ synced=true;
+ }
+};
+
+#else
/* Container for Body External Variables (bex), typically forces and torques from interactions.
* Values should be reset at every iteration by BexResetter.
* If you want to add your own bex, there are 4 steps:
@@ -22,21 +99,22 @@
std::vector<Vector3r> _force;
std::vector<Vector3r> _torque;
size_t size;
+ inline void ensureSize(body_id_t id){ if(size<=(size_t)id) resize(min((size_t)1.5*(id+100),(size_t)(id+2000)));}
public:
BexContainer(): size(0){}
- inline void ensureSize(body_id_t id){ if(size<=(size_t)id) resize(min((size_t)1.5*(id+100),(size_t)(id+2000)));}
- //! Get writable reference to force acting on body # id
- inline Vector3r& force(body_id_t id){ ensureSize(id); return _force[id]; }
- //! Get writable reference to torque acting on body # id
- inline Vector3r& torque(body_id_t id){ ensureSize(id); return _torque[id]; }
+ const Vector3r& getForce(body_id_t id){ensureSize(id); return _force[id];}
+ void addForce(body_id_t id,const Vector3r& f){ensureSize(id); _force[id]+=f;}
+ const Vector3r& getTorque(body_id_t id){ensureSize(id); return _torque[id];}
+ void addTorque(body_id_t id,const Vector3r& t){ensureSize(id); _torque[id]+=t;}
//! Set all bex's to zero
void reset(){
memset(_force[0], 0,sizeof(Vector3r)*size);
memset(_torque[0],0,sizeof(Vector3r)*size);
}
+ //! No-op for API compatibility with the threaded version
+ void sync(){return;}
/*! Resize the container; this happens automatically,
- * but you may want to set the size beforehand to avoid resizes as the simulation grows.
- */
+ * but you may want to set the size beforehand to avoid resizes as the simulation grows. */
void resize(size_t newSize){
_force.resize(newSize);
_torque.resize(newSize);
@@ -44,3 +122,6 @@
std::cerr<<"[DEBUG] BexContainer: Resized to "<<size<<std::endl;
}
};
+
+
+#endif
Modified: trunk/core/Omega.hpp
===================================================================
--- trunk/core/Omega.hpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/core/Omega.hpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -167,6 +167,7 @@
Omega& operator=(const Omega&);
FRIEND_SINGLETON(Omega);
+ friend class pyOmega;
};
Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/extra/Brefcom.cpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -12,7 +12,8 @@
CREATE_LOGGER(BrefcomGlobalCharacteristics);
void BrefcomGlobalCharacteristics::compute(MetaBody* rb, bool useMaxForce){
- Shop::Bex::initCache();
+ //Shop::Bex::initCache();
+ rb->bex.sync();
// 1. reset volumetric strain (cummulative in the next loop)
// 2. get maximum force on a body and sum of all forces (for averaging)
@@ -21,7 +22,7 @@
BrefcomPhysParams* bpp(YADE_CAST<BrefcomPhysParams*>(b->physicalParameters.get()));
bpp->epsVolumetric=0;
bpp->numContacts=0;
- currF=Shop::Bex::force(b->id,rb).Length(); maxF=max(currF,maxF); sumF+=currF;
+ currF=rb->bex.getForce(b->id).Length(); maxF=max(currF,maxF); sumF+=currF;
}
Real meanF=sumF/rb->bodies->size();
@@ -124,10 +125,10 @@
CREATE_LOGGER(BrefcomLaw);
void BrefcomLaw::applyForce(const Vector3r& force, const body_id_t& id1, const body_id_t& id2){
- Shop::Bex::force(id1,rootBody)+=force;
- Shop::Bex::force(id2,rootBody)-=force;
- Shop::Bex::momentum(id1,rootBody)+=(contGeom->contactPoint-contGeom->pos1).Cross(force);
- Shop::Bex::momentum(id2,rootBody)+=(contGeom->contactPoint-contGeom->pos2).Cross(-force);
+ rootBody->bex.addForce(id1,force);
+ rootBody->bex.addForce(id2,-force);
+ rootBody->bex.addTorque(id1,(contGeom->contactPoint-contGeom->pos1).Cross(force));
+ rootBody->bex.addTorque(id2,(contGeom->contactPoint-contGeom->pos2).Cross(-force));
}
CREATE_LOGGER(ef2_Spheres_Brefcom_BrefcomLaw);
Modified: trunk/extra/clump/Shop.cpp
===================================================================
--- trunk/extra/clump/Shop.cpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/extra/clump/Shop.cpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -39,7 +39,6 @@
#include<yade/pkg-common/Force.hpp>
#include<yade/pkg-common/Momentum.hpp>
-#include<yade/pkg-dem/GlobalStiffness.hpp>
/*class InteractingSphere2AABB;
class InteractingBox2AABB;
class MetaInteractingGeometry;
@@ -87,18 +86,16 @@
int Shop::Bex::forceIdx=-1;
int Shop::Bex::momentumIdx=-1;
-int Shop::Bex::globalStiffnessIdx=-1;
void Shop::Bex::initCache(){
if(Shop::Bex::forceIdx<0){
Shop::Bex::forceIdx=shared_ptr<PhysicalAction>(new Force())->getClassIndex();
Shop::Bex::momentumIdx=shared_ptr<PhysicalAction>(new Momentum())->getClassIndex();
- Shop::Bex::globalStiffnessIdx=shared_ptr<PhysicalAction>(new GlobalStiffness())->getClassIndex();
}
}
#ifdef BEX_CONTAINER
- Vector3r& Shop::Bex::force(body_id_t id,MetaBody* rb){ return rb->bex.force(id);}
- Vector3r& Shop::Bex::momentum(body_id_t id,MetaBody* rb){ return rb->bex.force(id);}
+ const Vector3r& Shop::Bex::force(body_id_t id,MetaBody* rb){ rb->bex.sync(); return rb->bex.getForce(id);}
+ const Vector3r& Shop::Bex::momentum(body_id_t id,MetaBody* rb){ rb->bex.sync(); return rb->bex.getTorque(id);}
#else
#define __BEX_ACCESS(retType,shopBexMember,bexClass,bexIdx,bexAttribute) retType& Shop::Bex::shopBexMember(body_id_t id,MetaBody* mb){ assert(bexIdx>=0); shared_ptr<PhysicalActionContainer> pac=(mb?mb:Omega::instance().getRootBody().get())->physicalActions; /*if((long)pac->size()<=id) throw invalid_argument("No " #shopBexMember " for #"+lexical_cast<string>(id)+" (max="+lexical_cast<string>(((long)pac->size()-1))+")");*/ return static_pointer_cast<bexClass>(pac->find(id,bexIdx))->bexAttribute; }
__BEX_ACCESS(Vector3r,force,Force,forceIdx,force);
@@ -110,10 +107,17 @@
*
* Shop::Bex::initCache must have been called at some point. */
void Shop::applyForceAtContactPoint(const Vector3r& force, const Vector3r& contPt, body_id_t id1, const Vector3r& pos1, body_id_t id2, const Vector3r& pos2, MetaBody* rootBody){
+#ifdef BEX_CONTAINER
+ rootBody->bex.addForce(id1,force);
+ rootBody->bex.addForce(id2,-force);
+ rootBody->bex.addTorque(id1,(contPt-pos1).Cross(force));
+ rootBody->bex.addTorque(id2,-(contPt-pos2).Cross(force));
+#else
Shop::Bex::force(id1,rootBody)+=force;
Shop::Bex::force(id2,rootBody)-=force;
Shop::Bex::momentum(id1,rootBody)+=(contPt-pos1).Cross(force);
Shop::Bex::momentum(id2,rootBody)-=(contPt-pos2).Cross(force);
+#endif
}
@@ -124,7 +128,12 @@
Real sumF=0,maxF=0,currF;
FOREACH(const shared_ptr<Body>& b, *rb->bodies){
if(!b->isDynamic) continue;
- currF=Shop::Bex::force(b->id,rb).Length(); maxF=max(currF,maxF); sumF+=currF;
+ #ifdef BEX_CONTAINER
+ LOG_FATAL("No implemented with BEX_CONTAINER yet");
+ abort();
+ #else
+ currF=Shop::Bex::force(b->id,rb).Length(); maxF=max(currF,maxF); sumF+=currF;
+ #endif
}
Real meanF=sumF/rb->bodies->size();
// get max force on contacts
Modified: trunk/extra/clump/Shop.hpp
===================================================================
--- trunk/extra/clump/Shop.hpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/extra/clump/Shop.hpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -86,12 +86,15 @@
*/
class Bex{
public:
- static int forceIdx,momentumIdx,globalStiffnessIdx;
+ static int forceIdx,momentumIdx;
static void initCache();
- static Vector3r& force(body_id_t, MetaBody* mb=NULL);
- static Vector3r& momentum(body_id_t, MetaBody* mb=NULL);
- static Vector3r& globalStiffness(body_id_t, MetaBody* mb=NULL);
- static Vector3r& globalRStiffness(body_id_t, MetaBody* mb=NULL);
+ #ifdef BEX_CONTAINER
+ static const Vector3r& force(body_id_t, MetaBody* mb=NULL);
+ static const Vector3r& momentum(body_id_t, MetaBody* mb=NULL);
+ #else
+ static Vector3r& force(body_id_t, MetaBody* mb=NULL);
+ static Vector3r& momentum(body_id_t, MetaBody* mb=NULL);
+ #endif
};
//! Estimate timestep based on P-wave propagation speed
Modified: trunk/extra/usct/UniaxialStrainControlledTest.cpp
===================================================================
--- trunk/extra/usct/UniaxialStrainControlledTest.cpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/extra/usct/UniaxialStrainControlledTest.cpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -129,8 +129,9 @@
void UniaxialStrainer::computeAxialForce(MetaBody* rootBody){
sumPosForces=sumNegForces=0;
#ifdef BEX_CONTAINER
- FOREACH(body_id_t id, negIds) sumNegForces+=rootBody->bex.force(id)[axis];
- FOREACH(body_id_t id, posIds) sumNegForces-=rootBody->bex.force(id)[axis];
+ rootBody->bex.sync();
+ FOREACH(body_id_t id, negIds) sumNegForces+=rootBody->bex.getForce(id)[axis];
+ FOREACH(body_id_t id, posIds) sumNegForces-=rootBody->bex.getForce(id)[axis];
#else
FOREACH(body_id_t id, negIds) sumNegForces+=Shop::Bex::force(id)[axis];
FOREACH(body_id_t id, posIds) sumPosForces-=Shop::Bex::force(id)[axis];
Modified: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/gui/py/_utils.cpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -187,14 +187,19 @@
* projected onto axis.
*/
Real sumBexTorques(int mask, python::tuple _axis, python::tuple _axisPt){
+ shared_ptr<MetaBody> rb=Omega::instance().getRootBody();
Shop::Bex::initCache();
- shared_ptr<MetaBody> rb=Omega::instance().getRootBody();
- Real ret;
+ Real ret=0;
Vector3r axis=tuple2vec(_axis), axisPt=tuple2vec(_axisPt);
FOREACH(const shared_ptr<Body> b, *rb->bodies){
if(!b->maskOk(mask)) continue;
- const Vector3r& m=Shop::Bex::momentum(b->getId(),rb.get());
- const Vector3r& f=Shop::Bex::force(b->getId(),rb.get());
+ #ifdef BEX_CONTAINER
+ const Vector3r& m=rb->bex.getTorque(b->getId());
+ const Vector3r& f=rb->bex.getForce(b->getId());
+ #else
+ const Vector3r& m=Shop::Bex::momentum(b->getId(),rb.get());
+ const Vector3r& f=Shop::Bex::force(b->getId(),rb.get());
+ #endif
Vector3r r=b->physicalParameters->se3.position-axisPt;
ret+=axis.Dot(m+r.Cross(f));
}
@@ -213,7 +218,11 @@
Vector3r direction=tuple2vec(_direction);
FOREACH(const shared_ptr<Body> b, *rb->bodies){
if(!b->maskOk(mask)) continue;
- const Vector3r& f=Shop::Bex::force(b->getId(),rb.get());
+ #ifdef BEX_CONTAINER
+ const Vector3r& f=rb->bex.getForce(b->getId());
+ #else
+ const Vector3r& f=Shop::Bex::force(b->getId(),rb.get());
+ #endif
ret+=direction.Dot(f);
}
return ret;
Modified: trunk/gui/py/yadeControl.cpp
===================================================================
--- trunk/gui/py/yadeControl.cpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/gui/py/yadeControl.cpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -387,8 +387,8 @@
class pyBexContainer{
public:
pyBexContainer(){}
- python::tuple force_get(long id){ Vector3r& f=Omega::instance().getRootBody()->bex.force(id); return python::make_tuple(f[0],f[1],f[2]); }
- python::tuple torque_get(long id){ Vector3r& m=Omega::instance().getRootBody()->bex.torque(id); return python::make_tuple(m[0],m[1],m[2]);}
+ python::tuple force_get(long id){ MetaBody* rb=Omega::instance().getRootBody().get(); rb->bex.sync(); Vector3r f=rb->bex.getForce(id); return python::make_tuple(f[0],f[1],f[2]); }
+ python::tuple torque_get(long id){ MetaBody* rb=Omega::instance().getRootBody().get(); rb->bex.sync(); Vector3r m=rb->bex.getTorque(id); return python::make_tuple(m[0],m[1],m[2]);}
};
#endif
@@ -459,13 +459,13 @@
void loadTmp(string mark){ load(":memory:"+mark);}
void tmpToFile(string mark, string filename){
// FIXME: memSavedSimulations are private, but I don't want to recompile all yade now; move it to public and uncomment these few lines at some point
- // if(OMEGA.memSavedSimulations.count(":memory:"+mark)==0) throw runtime_error("No memory-saved simulation named "+mark);
+ if(OMEGA.memSavedSimulations.count(":memory:"+mark)==0) throw runtime_error("No memory-saved simulation named "+mark);
iostreams::filtering_ostream out;
if(boost::algorithm::ends_with(filename,".bz2")) out.push(iostreams::bzip2_compressor());
out.push(iostreams::file_sink(filename));
if(!out.good()) throw runtime_error("Error while opening file `"+filename+"' for writing.");
LOG_INFO("Saving :memory:"<<mark<<" to "<<filename);
- //out<<OMEGA.memSavedSimulations[":memory:"+mark];
+ out<<OMEGA.memSavedSimulations[":memory:"+mark];
}
Modified: trunk/pkg/common/Engine/DeusExMachina/GravityEngines.cpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/GravityEngines.cpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/pkg/common/Engine/DeusExMachina/GravityEngines.cpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -25,7 +25,7 @@
shared_ptr<ParticleParameters> p=YADE_PTR_CAST<ParticleParameters>(b->physicalParameters);
if(p!=0) //not everything derives from ParticleParameters; this line was assert(p); - Janek
#ifdef BEX_CONTAINER
- ncb->bex.force(b->getId())+=gravity*p->mass;
+ ncb->bex.addForce(b->getId(),gravity*p->mass);
#else
static_cast<Force*>(ncb->physicalActions->find(b->getId(),cachedForceClassIndex).get())->force+=gravity*p->mass;
#endif
@@ -39,8 +39,8 @@
Real F=accel*YADE_PTR_CAST<ParticleParameters>(b->physicalParameters)->mass;
Vector3r toCenter=centralPos-b->physicalParameters->se3.position; toCenter.Normalize();
#ifdef BEX_CONTAINER
- rootBody->bex.force(b->getId())+=F*toCenter;
- if(reciprocal) rootBody->bex.force(centralBody)-=F*toCenter;
+ rootBody->bex.addForce(b->getId(),F*toCenter);
+ if(reciprocal) rootBody->bex.addForce(centralBody,-F*toCenter);
#else
static_cast<Force*>(rootBody->physicalActions->find(b->getId(),cachedForceClassIndex).get())->force+=F*toCenter;
if(reciprocal) static_cast<Force*>(rootBody->physicalActions->find(centralBody,cachedForceClassIndex).get())->force-=F*toCenter;
@@ -59,7 +59,7 @@
Vector3r closestAxisPoint=(x2-x1) * /* t */ (-(x1-x0).Dot(x2-x1))/((x2-x1).SquaredLength());
Vector3r toAxis=closestAxisPoint-x0; toAxis.Normalize();
#ifdef BEX_CONTAINER
- rootBody->bex.force(b->getId())+=acceleration*myMass*toAxis;
+ rootBody->bex.addForce(b->getId(),acceleration*myMass*toAxis);
#else
static_pointer_cast<Force>(rootBody->physicalActions->find(b->getId(),cachedForceClassIndex))->force+=acceleration*myMass*toAxis;
//if(b->getId()==20){ TRVAR3(toAxis,acceleration*myMass*toAxis,static_pointer_cast<Force>(rootBody->physicalActions->find(b->getId(),cachedForceClassIndex))->force); }
Modified: trunk/pkg/common/Engine/MetaEngine/ConstitutiveLaw.hpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/ConstitutiveLaw.hpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/pkg/common/Engine/MetaEngine/ConstitutiveLaw.hpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -23,8 +23,8 @@
* is done only if the derived ConstitutiveLaw says NEEDS_BEX("Force","Momentum"), for example.
*/
#ifdef BEX_CONTAINER
- inline Vector3r& bodyForce(const body_id_t id, MetaBody* rb) const { return rb->bex.force(id); }
- inline Vector3r& bodyTorque(const body_id_t id, MetaBody* rb) const { return rb->bex.torque(id);}
+ void addForce (const body_id_t id, const Vector3r& f,MetaBody* rb){rb->bex.addForce (id,f);}
+ void addTorque(const body_id_t id, const Vector3r& t,MetaBody* rb){rb->bex.addTorque(id,t);}
#else
Vector3r& bodyForce(const body_id_t id, MetaBody* rb) const { return static_pointer_cast<Force>(rb->physicalActions->find(id,forceIdx))->force; }
Vector3r& bodyTorque(const body_id_t id, MetaBody* rb) const { return static_pointer_cast<Momentum>(rb->physicalActions->find(id,torqueIdx))->momentum;}
@@ -33,9 +33,9 @@
/*! Convenience function to apply force and torque from one force at contact point. Not sure if this is the right place for it. */
void applyForceAtContactPoint(const Vector3r& force, const Vector3r& contactPoint, const body_id_t id1, const Vector3r& pos1, const body_id_t id2, const Vector3r& pos2, MetaBody* rb){
#ifdef BEX_CONTAINER
- rb->bex.force(id1)+=force; rb->bex.force(id2)-=force;
- rb->bex.torque(id1)+=(contactPoint-pos1).Cross(force);
- rb->bex.torque(id2)-=(contactPoint-pos2).Cross(force);
+ rb->bex.addForce(id1,force); rb->bex.addForce(id2,-force);
+ rb->bex.addTorque(id1,(contactPoint-pos1).Cross(force));
+ rb->bex.addTorque(id2,-(contactPoint-pos2).Cross(force));
#else
bodyForce(id1,rb)+=force; bodyForce(id2,rb)-=force;
bodyTorque(id1,rb)+=(contactPoint-pos1).Cross(force);
Modified: trunk/pkg/common/Engine/StandAloneEngine/PhysicalActionContainerInitializer.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/PhysicalActionContainerInitializer.cpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/pkg/common/Engine/StandAloneEngine/PhysicalActionContainerInitializer.cpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -32,6 +32,9 @@
void PhysicalActionContainerInitializer::action(MetaBody* ncb)
{
#ifdef BEX_CONTAINER
+ // this is not necessary, since BexContainer grows upon requests.
+ // But it takes about 10 resizes to get to 2000 bodies (only at the beginning), so you can save a few milliseconds here
+ ncb->bex.resize(ncb->bodies->size());
return;
#endif
list<string> allNames;
Modified: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -31,17 +31,22 @@
void NewtonsDampedLaw::applyCondition ( MetaBody * ncb )
{
+ #ifdef BEX_CONTAINER
+ #ifdef YADE_OPENMP
+ ncb->bex.sync();
+ #endif
+ #endif
FOREACH(const shared_ptr<Body>& b, *ncb->bodies){
if (!b->isDynamic) continue;
RigidBodyParameters* rb = YADE_CAST<RigidBodyParameters*>(b->physicalParameters.get());
unsigned int id = b->getId();
#ifdef BEX_CONTAINER
- Vector3r& m = ncb->bex.torque(id);
- Vector3r& f = ncb->bex.force(id);
+ const Vector3r& m=ncb->bex.getTorque(id);
+ const Vector3r& f=ncb->bex.getForce(id);
#else
- Vector3r& m = ( static_cast<Momentum*> ( ncb->physicalActions->find ( id, momentumClassIndex ).get() ) )->momentum;
- Vector3r& f = ( static_cast<Force*> ( ncb->physicalActions->find ( id, forceClassIndex ).get() ) )->force;
+ const Vector3r& m = ( static_cast<Momentum*> ( ncb->physicalActions->find ( id, momentumClassIndex ).get() ) )->momentum;
+ const Vector3r& f = ( static_cast<Force*> ( ncb->physicalActions->find ( id, forceClassIndex ).get() ) )->force;
#endif
Real dt = Omega::instance().getTimeStep();
Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -198,7 +198,10 @@
{
//cerr << "TriaxialStressController::applyCondition" << endl;
+ // sync thread storage of BexContainer
+ ncb->bex.sync();
+
shared_ptr<BodyContainer>& bodies = ncb->bodies;
Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.hpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -33,7 +33,7 @@
bool firstRun;
inline const Vector3r getForce(MetaBody* rb, body_id_t id){
#ifdef BEX_CONTAINER
- return rb->bex.force(id);
+ return rb->bex.getForce(id); // needs sync, which is done at the beginning of applyCondition
#else
return static_cast<Force*>(rb->physicalActions->find(id,ForceClassIndex).get())->force
#endif
Modified: trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -67,8 +67,8 @@
}
Vector3r f = phys->normalForce + shearForce;
- bodyForce (id1,rootBody) -= f;
- bodyForce (id2,rootBody) += f;
- bodyTorque(id1,rootBody) -= c1x.Cross(f);
- bodyTorque(id2,rootBody) += c2x.Cross(f);
+ addForce (id1,-f,rootBody);
+ addForce (id2, f,rootBody);
+ addTorque(id1,-c1x.Cross(f),rootBody);
+ addTorque(id2, c2x.Cross(f),rootBody);
}
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-02-22 13:08:47 UTC (rev 1679)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-02-22 20:08:45 UTC (rev 1680)
@@ -44,10 +44,15 @@
Shop::applyForceAtContactPoint(force,contGeom->contactPoint,i->getId1(),contGeom->pos1,i->getId2(),contGeom->pos2,rb);
Vector3r bendAbs; Real torsionAbs; contGeom->bendingTorsionAbs(bendAbs,torsionAbs);
- Shop::Bex::momentum(i->getId1(),rb)+=contGeom->normal*torsionAbs*ktor;
- Shop::Bex::momentum(i->getId2(),rb)-=contGeom->normal*torsionAbs*ktor;
- Shop::Bex::momentum(i->getId1(),rb)+=bendAbs*kb;
- Shop::Bex::momentum(i->getId2(),rb)-=bendAbs*kb;
+ #ifdef BEX_CONTAINER
+ rb->bex.addTorque(i->getId1(), contGeom->normal*torsionAbs*ktor+bendAbs*kb);
+ rb->bex.addTorque(i->getId2(),-contGeom->normal*torsionAbs*ktor-bendAbs*kb);
+ #else
+ Shop::Bex::momentum(i->getId1(),rb)+=contGeom->normal*torsionAbs*ktor;
+ Shop::Bex::momentum(i->getId2(),rb)-=contGeom->normal*torsionAbs*ktor;
+ Shop::Bex::momentum(i->getId1(),rb)+=bendAbs*kb;
+ Shop::Bex::momentum(i->getId2(),rb)-=bendAbs*kb;
+ #endif
}
}
@@ -170,10 +175,10 @@
Vector3r f = currentContactPhysics->normalForce + shearForce;
#ifdef BEX_CONTAINER
- ncb->bex.force(id1)-=f;
- ncb->bex.force(id2)+=f;
- ncb->bex.torque(id1)-=c1x.Cross(f);
- ncb->bex.torque(id2)+=c2x.Cross(f);
+ ncb->bex.addForce (id1,-f);
+ ncb->bex.addForce (id2,+f);
+ ncb->bex.addTorque(id1,-c1x.Cross(f));
+ ncb->bex.addTorque(id2, c2x.Cross(f));
#else
static_cast<Force*> ( ncb->physicalActions->find( id1 , actionForceIndex).get() )->force -= f;
static_cast<Force*> ( ncb->physicalActions->find( id2 , actionForceIndex ).get() )->force += f;
