yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01465
[svn] r1863 - in trunk: core core/containers pkg/common pkg/common/DataClass pkg/common/Engine/MetaEngine pkg/dem/Engine/DeusExMachina scripts
Author: eudoxos
Date: 2009-07-13 23:39:27 +0200 (Mon, 13 Jul 2009)
New Revision: 1863
Modified:
trunk/core/BodyContainer.hpp
trunk/core/containers/BodyAssocVector.cpp
trunk/core/containers/BodyAssocVector.hpp
trunk/core/containers/BodyRedirectionVector.cpp
trunk/core/containers/BodyRedirectionVector.hpp
trunk/pkg/common/DataClass/VelocityBins.cpp
trunk/pkg/common/DataClass/VelocityBins.hpp
trunk/pkg/common/Engine/MetaEngine/EngineUnits.cpp
trunk/pkg/common/Engine/MetaEngine/PhysicalActionApplierUnit.hpp
trunk/pkg/common/Engine/MetaEngine/PhysicalActionDamperUnit.hpp
trunk/pkg/common/SConscript
trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp
trunk/scripts/simple-scene.py
Log:
1. Parallelized NewtonsDampedLaw (scales very well)
2. Parallelized VelocityBins
3. Add const modifier to BodyContainer::size
4. Add virtual dtors to PhysicalAction{Damper,Applier}Unit (for dynamic_cast'ing)
Modified: trunk/core/BodyContainer.hpp
===================================================================
--- trunk/core/BodyContainer.hpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/core/BodyContainer.hpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -87,14 +87,14 @@
virtual bool erase(unsigned int) {throw;};
virtual bool find(unsigned int , shared_ptr<Body>&) const {throw;};
virtual bool exists(unsigned int id) const {throw;};
- virtual shared_ptr<Body>& operator[](unsigned int) {throw;};
- virtual const shared_ptr<Body>& operator[](unsigned int) const {throw;};
+ virtual shared_ptr<Body>& operator[](unsigned int) {throw;};
+ virtual const shared_ptr<Body>& operator[](unsigned int) const {throw;};
typedef BodyContainerIteratorPointer iterator;
virtual BodyContainer::iterator begin() {throw;};
virtual BodyContainer::iterator end() {throw;};
- virtual unsigned int size() {throw;};
+ virtual unsigned int size() const {throw;};
protected :
void setId(shared_ptr<Body>& , unsigned int);
Modified: trunk/core/containers/BodyAssocVector.cpp
===================================================================
--- trunk/core/containers/BodyAssocVector.cpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/core/containers/BodyAssocVector.cpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -210,7 +210,7 @@
// }
//
-unsigned int BodyAssocVector::size()
+unsigned int BodyAssocVector::size() const
{
return bodies.size();
}
Modified: trunk/core/containers/BodyAssocVector.hpp
===================================================================
--- trunk/core/containers/BodyAssocVector.hpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/core/containers/BodyAssocVector.hpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -60,7 +60,7 @@
virtual BodyContainer::iterator begin();
virtual BodyContainer::iterator end();
- virtual unsigned int size();
+ virtual unsigned int size() const;
REGISTER_CLASS_NAME(BodyAssocVector);
REGISTER_BASE_CLASS_NAME(BodyContainer);
Modified: trunk/core/containers/BodyRedirectionVector.cpp
===================================================================
--- trunk/core/containers/BodyRedirectionVector.cpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/core/containers/BodyRedirectionVector.cpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -255,7 +255,7 @@
// }
//
-unsigned int BodyRedirectionVector::size()
+unsigned int BodyRedirectionVector::size() const
{
return bodies.size();
}
Modified: trunk/core/containers/BodyRedirectionVector.hpp
===================================================================
--- trunk/core/containers/BodyRedirectionVector.hpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/core/containers/BodyRedirectionVector.hpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -55,7 +55,7 @@
virtual BodyContainer::iterator begin();
virtual BodyContainer::iterator end();
- virtual unsigned int size();
+ virtual unsigned int size() const;
// serialization of this class...
REGISTER_CLASS_NAME(BodyRedirectionVector);
Modified: trunk/pkg/common/DataClass/VelocityBins.cpp
===================================================================
--- trunk/pkg/common/DataClass/VelocityBins.cpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/pkg/common/DataClass/VelocityBins.cpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -38,7 +38,7 @@
if(refMaxVelSq<0){ refMaxVelSq=currMaxVelSq; /* first time */}
else {
// there should be some maximum speed change parameter, so that bins do not change their limits (and therefore bodies, also!) too often, depending on 1 particle going crazy
- refMaxVelSq=min(max(refMaxVelSq*pow(1-maxRefRelStep,2),currMaxVelSq),refMaxVelSq*pow(1+maxRefRelStep,2));
+ refMaxVelSq=min(max(refMaxVelSq/pow(1+maxRefRelStep,2),currMaxVelSq),refMaxVelSq*pow(1+maxRefRelStep,2));
if(refMaxVelSq==0) refMaxVelSq=currMaxVelSq;
}
LOG_TRACE("new refMaxVel: "<<sqrt(refMaxVelSq));
@@ -86,12 +86,23 @@
}
/* non-parallel implementations */
-void VelocityBins::binVelSqInitialize(){ FOREACH(Bin& bin, bins) bin.currMaxVelSq=0; }
-void VelocityBins::binVelSqUse(body_id_t id, Real velSq){
- #ifdef YADE_OPENMP
- if(omp_get_thread_num()>0) {LOG_FATAL("VelocityBins do not support multiple openMP threads yet!"); abort(); }
- #endif
- Real& maxVelSq(bins[bodyBins[id]].currMaxVelSq);
- maxVelSq=max(maxVelSq,velSq);
-}
-void VelocityBins::binVelSqFinalize(){}
+#ifdef YADE_OPENMP
+ void VelocityBins::binVelSqInitialize(){ FOREACH(Bin& bin, bins){ FOREACH(Real& vSq, bin.threadMaxVelSq) vSq=0.; }}
+ void VelocityBins::binVelSqUse(body_id_t id, Real velSq){
+ Real& maxVelSq(bins[bodyBins[id]].threadMaxVelSq[omp_get_thread_num()]);
+ maxVelSq=max(maxVelSq,velSq);
+ }
+ void VelocityBins::binVelSqFinalize(){
+ FOREACH(Bin& bin, bins){
+ bin.currMaxVelSq=0;
+ FOREACH(const Real& vSq, bin.threadMaxVelSq) bin.currMaxVelSq=max(bin.currMaxVelSq,vSq);
+ }
+ }
+#else
+ void VelocityBins::binVelSqInitialize(){ FOREACH(Bin& bin, bins) bin.currMaxVelSq=0; }
+ void VelocityBins::binVelSqUse(body_id_t id, Real velSq){
+ Real& maxVelSq(bins[bodyBins[id]].currMaxVelSq);
+ maxVelSq=max(maxVelSq,velSq);
+ }
+ void VelocityBins::binVelSqFinalize(){}
+#endif
Modified: trunk/pkg/common/DataClass/VelocityBins.hpp
===================================================================
--- trunk/pkg/common/DataClass/VelocityBins.hpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/pkg/common/DataClass/VelocityBins.hpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -4,6 +4,9 @@
#include<yade/core/Interaction.hpp> // for body_id_t
#include<yade/pkg-common/RigidBodyParameters.hpp>
#include<vector>
+#ifdef YADE_OPENMP
+ #include<omp.h>
+#endif
class MetaBody;
@@ -14,10 +17,14 @@
*/
class VelocityBins{
public:
- VelocityBins(int _nBins, Real _refMaxVelSq, Real _binCoeff=10, Real _binOverlap=0.8): refMaxVelSq(_refMaxVelSq), binCoeff(_binCoeff), binOverlap(_binOverlap), maxRefRelStep(0.05), nBins(_nBins), histInterval(200), histLast(-1){};
+ VelocityBins(int _nBins, Real _refMaxVelSq, Real _binCoeff=10, Real _binOverlap=0.8): refMaxVelSq(_refMaxVelSq), binCoeff(_binCoeff), binOverlap(_binOverlap), maxRefRelStep(100), nBins(_nBins), histInterval(200), histLast(-1){}
typedef signed char binNo_t;
struct Bin{
- Bin(): binMinVelSq(-1), binMaxVelSq(-1), maxDist(0), currDistSq(0), currMaxVelSq(0), nBodies(0){}
+ Bin(): binMinVelSq(-1), binMaxVelSq(-1), maxDist(0), currDistSq(0), currMaxVelSq(0), nBodies(0){
+ #ifdef YADE_OPENMP
+ threadMaxVelSq.resize(omp_get_max_threads());
+ #endif
+ };
// limits for bin memebrship
Real binMinVelSq, binMaxVelSq;
// maximum distance that body in this bin can travel before it goes out of its swept bbox
@@ -28,6 +35,9 @@
Real currMaxVelSq;
// number of bodies in this bin (for informational purposes only)
long nBodies;
+ #ifdef YADE_OPENMP
+ vector<Real> threadMaxVelSq;
+ #endif
};
// bins themselves (their number is nBins)
vector<Bin> bins;
Modified: trunk/pkg/common/Engine/MetaEngine/EngineUnits.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/EngineUnits.cpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/pkg/common/Engine/MetaEngine/EngineUnits.cpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -4,6 +4,8 @@
#include<yade/pkg-common/InteractionGeometryEngineUnit.hpp>
#include<yade/pkg-common/InteractionPhysicsEngineUnit.hpp>
#include<yade/pkg-common/PhysicalParametersEngineUnit.hpp>
+#include<yade/pkg-common/PhysicalActionDamperUnit.hpp>
+#include<yade/pkg-common/PhysicalActionApplierUnit.hpp>
#include<yade/pkg-common/ConstitutiveLaw.hpp>
BoundingVolumeEngineUnit::~BoundingVolumeEngineUnit(){};
@@ -12,5 +14,7 @@
InteractionGeometryEngineUnit::~InteractionGeometryEngineUnit(){};
InteractionPhysicsEngineUnit::~InteractionPhysicsEngineUnit(){};
PhysicalParametersEngineUnit::~PhysicalParametersEngineUnit(){};
+PhysicalActionDamperUnit::~PhysicalActionDamperUnit(){};
+PhysicalActionApplierUnit::~PhysicalActionApplierUnit(){};
ConstitutiveLaw::~ConstitutiveLaw(){};
Modified: trunk/pkg/common/Engine/MetaEngine/PhysicalActionApplierUnit.hpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/PhysicalActionApplierUnit.hpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/pkg/common/Engine/MetaEngine/PhysicalActionApplierUnit.hpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -14,6 +14,7 @@
#include<yade/core/EngineUnit1D.hpp>
class PhysicalActionApplierUnit: public EngineUnit1D<void,TYPELIST_3(const shared_ptr<PhysicalParameters>&,const Body*, MetaBody*)>{
+ public: virtual ~PhysicalActionApplierUnit();
REGISTER_CLASS_AND_BASE(PhysicalActionApplierUnit,EngineUnit1D);
REGISTER_ATTRIBUTES(EngineUnit, /* nothing here */ );
};
Modified: trunk/pkg/common/Engine/MetaEngine/PhysicalActionDamperUnit.hpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/PhysicalActionDamperUnit.hpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/pkg/common/Engine/MetaEngine/PhysicalActionDamperUnit.hpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -17,6 +17,8 @@
class PhysicalActionDamperUnit: public EngineUnit1D<void,TYPELIST_3(const shared_ptr<PhysicalParameters>&,const Body*, MetaBody*)>{
REGISTER_CLASS_AND_BASE(PhysicalActionDamperUnit,EngineUnit1D);
REGISTER_ATTRIBUTES(EngineUnit,/* nothing here */);
+ public: virtual ~PhysicalActionDamperUnit();
+ protected:
/* We are friend of BexContainer. These functions can be used safely provided that bex is NEVER read after being modified. */
Vector3r getForceUnsynced (body_id_t id, MetaBody* rb){ return rb->bex.getForceUnsynced (id);}
Vector3r getTorqueUnsynced(body_id_t id, MetaBody* rb){ return rb->bex.getTorqueUnsynced(id);}
Modified: trunk/pkg/common/SConscript
===================================================================
--- trunk/pkg/common/SConscript 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/pkg/common/SConscript 2009-07-13 21:39:27 UTC (rev 1863)
@@ -97,7 +97,7 @@
LIBS=env['LIBS']+['FilterEngine','ColorScale']),
env.SharedLibrary('ColorizedTimeFilter',['Engine/FilterEngine/ColorizedTimeFilter.cpp'],
LIBS=env['LIBS']+['FilterEngine','ColorScale']),
- env.SharedLibrary('BoundingVolumeMetaEngine',['Engine/MetaEngine/BoundingVolumeMetaEngine.cpp'],LIBS=env['LIBS']+['RigidBodyParameters','AABB','EngineUnits']),
+ env.SharedLibrary('BoundingVolumeMetaEngine',['Engine/MetaEngine/BoundingVolumeMetaEngine.cpp'],LIBS=env['LIBS']+['RigidBodyParameters','AABB','EngineUnits','VelocityBins']),
env.SharedLibrary('GeometricalModelMetaEngine',['Engine/MetaEngine/GeometricalModelMetaEngine.cpp'],LIBS=env['LIBS']+['EngineUnits']),
env.SharedLibrary('InteractingGeometryMetaEngine',['Engine/MetaEngine/InteractingGeometryMetaEngine.cpp'],LIBS=env['LIBS']+['EngineUnits']),
env.SharedLibrary('PhysicalParametersMetaEngine',['Engine/MetaEngine/PhysicalParametersMetaEngine.cpp'],LIBS=env['LIBS']+['EngineUnits']),
Modified: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -36,7 +36,11 @@
clumpRBP->acceleration+=diffClumpAccel;
clumpRBP->angularAcceleration+=diffClumpAngularAccel;
if(haveBins) velocityBins->binVelSqUse(memberId,VelocityBins::getBodyVelSq(rb));
- maxVelocitySq=max(maxVelocitySq,rb->velocity.SquaredLength());
+ #ifdef YADE_OPENMP
+ Real& thrMaxVSq=threadMaxVelocitySq[omp_get_thread_num()]; thrMaxVSq=max(thrMaxVSq,rb->velocity.SquaredLength());
+ #else
+ maxVelocitySq=max(maxVelocitySq,rb->velocity.SquaredLength());
+ #endif
}
void NewtonsDampedLaw::action(MetaBody * ncb)
@@ -47,59 +51,85 @@
haveBins=(bool)velocityBins;
if(haveBins) velocityBins->binVelSqInitialize();
- FOREACH(const shared_ptr<Body>& b, *ncb->bodies){
- // clump members are non-dynamic; they skip the rest of loop once their forces are properly taken into account, however
- if (!b->isDynamic || b->isClumpMember()) continue;
- RigidBodyParameters* rb = YADE_CAST<RigidBodyParameters*>(b->physicalParameters.get());
- const body_id_t& id=b->getId();
- const Vector3r& m=ncb->bex.getTorque(id); const Vector3r& f=ncb->bex.getForce(id);
+ #ifdef YADE_OPENMP
+ FOREACH(Real& thrMaxVSq, threadMaxVelocitySq) { thrMaxVSq=0; }
+ const BodyContainer& bodies=*(ncb->bodies.get());
+ const long size=(long)bodies.size();
+ #pragma omp parallel for schedule(static)
+ for(long _id=0; _id<size; _id++){
+ const shared_ptr<Body>& b(bodies[_id]);
+ #else
+ FOREACH(const shared_ptr<Body>& b, *ncb->bodies){
+ #endif
+ RigidBodyParameters* rb = YADE_CAST<RigidBodyParameters*>(b->physicalParameters.get());
+ const body_id_t& id=b->getId();
+ // clump members are non-dynamic; we only get their velocities here
+ if (!b->isDynamic || b->isClumpMember()){
+ // FIXME: duplicated code from below; awaits https://bugs.launchpad.net/yade/+bug/398089 to be solved
+ if(haveBins) {velocityBins->binVelSqUse(id,VelocityBins::getBodyVelSq(rb));}
+ #ifdef YADE_OPENMP
+ Real& thrMaxVSq=threadMaxVelocitySq[omp_get_thread_num()]; thrMaxVSq=max(thrMaxVSq,rb->velocity.SquaredLength());
+ #else
+ maxVelocitySq=max(maxVelocitySq,rb->velocity.SquaredLength());
+ #endif
+ continue;
+ }
+ const Vector3r& m=ncb->bex.getTorque(id); const Vector3r& f=ncb->bex.getForce(id);
- if (b->isStandalone()){
- rb->acceleration=f/rb->mass;
- rb->angularAcceleration=diagDiv(m,rb->inertia);
- cundallDamp(dt,f,rb->velocity,rb->acceleration,m,rb->angularVelocity,rb->angularAcceleration);
- }
- else if (b->isClump()){
- rb->acceleration=rb->angularAcceleration=Vector3r::ZERO; // to make sure; should be reset in Clump::moveMembers
- FOREACH(Clump::memberMap::value_type mm, static_cast<Clump*>(b.get())->members){
- handleClumpMember(ncb,mm.first,rb);
+ if (b->isStandalone()){
+ rb->acceleration=f/rb->mass;
+ rb->angularAcceleration=diagDiv(m,rb->inertia);
+ cundallDamp(dt,f,rb->velocity,rb->acceleration,m,rb->angularVelocity,rb->angularAcceleration);
}
- // at this point, forces from clump members are already summed up, this is just for forces applied to clump proper, if there are such
- Vector3r dLinAccel=f/rb->mass, dAngAccel=diagDiv(m,rb->inertia);
- cundallDamp(dt,f,rb->velocity,dLinAccel,m,rb->angularVelocity,dAngAccel);
- rb->acceleration+=dLinAccel;
- rb->angularAcceleration+=dAngAccel;
- }
+ else if (b->isClump()){
+ rb->acceleration=rb->angularAcceleration=Vector3r::ZERO; // to make sure; should be reset in Clump::moveMembers
+ // sum force on clump memebrs, add them to the clump itself
+ FOREACH(Clump::memberMap::value_type mm, static_cast<Clump*>(b.get())->members){
+ handleClumpMember(ncb,mm.first,rb);
+ }
+ // at this point, forces from clump members are already summed up, this is just for forces applied to clump proper, if there are such
+ Vector3r dLinAccel=f/rb->mass, dAngAccel=diagDiv(m,rb->inertia);
+ cundallDamp(dt,f,rb->velocity,dLinAccel,m,rb->angularVelocity,dAngAccel);
+ rb->acceleration+=dLinAccel;
+ rb->angularAcceleration+=dAngAccel;
+ }
- if(haveBins) {velocityBins->binVelSqUse(id,VelocityBins::getBodyVelSq(rb));}
- maxVelocitySq=max(maxVelocitySq,rb->velocity.SquaredLength());
+ // blocking DOFs
+ if(rb->blockedDOFs==0){ /* same as: rb->blockedDOFs==PhysicalParameters::DOF_NONE */
+ rb->angularVelocity=rb->angularVelocity+dt*rb->angularAcceleration;
+ rb->velocity=rb->velocity+dt*rb->acceleration;
+ } else {
+ if((rb->blockedDOFs & PhysicalParameters::DOF_X)==0) rb->velocity[0]+=dt*rb->acceleration[0];
+ if((rb->blockedDOFs & PhysicalParameters::DOF_Y)==0) rb->velocity[1]+=dt*rb->acceleration[1];
+ if((rb->blockedDOFs & PhysicalParameters::DOF_Z)==0) rb->velocity[2]+=dt*rb->acceleration[2];
+ if((rb->blockedDOFs & PhysicalParameters::DOF_RX)==0) rb->angularVelocity[0]+=dt*rb->angularAcceleration[0];
+ if((rb->blockedDOFs & PhysicalParameters::DOF_RY)==0) rb->angularVelocity[1]+=dt*rb->angularAcceleration[1];
+ if((rb->blockedDOFs & PhysicalParameters::DOF_RZ)==0) rb->angularVelocity[2]+=dt*rb->angularAcceleration[2];
+ }
- // blocking DOFs
- if(rb->blockedDOFs==0){ /* same as: rb->blockedDOFs==PhysicalParameters::DOF_NONE */
- rb->angularVelocity=rb->angularVelocity+dt*rb->angularAcceleration;
- rb->velocity=rb->velocity+dt*rb->acceleration;
- } else {
- if((rb->blockedDOFs & PhysicalParameters::DOF_X)==0) rb->velocity[0]+=dt*rb->acceleration[0];
- if((rb->blockedDOFs & PhysicalParameters::DOF_Y)==0) rb->velocity[1]+=dt*rb->acceleration[1];
- if((rb->blockedDOFs & PhysicalParameters::DOF_Z)==0) rb->velocity[2]+=dt*rb->acceleration[2];
- if((rb->blockedDOFs & PhysicalParameters::DOF_RX)==0) rb->angularVelocity[0]+=dt*rb->angularAcceleration[0];
- if((rb->blockedDOFs & PhysicalParameters::DOF_RY)==0) rb->angularVelocity[1]+=dt*rb->angularAcceleration[1];
- if((rb->blockedDOFs & PhysicalParameters::DOF_RZ)==0) rb->angularVelocity[2]+=dt*rb->angularAcceleration[2];
- }
+ // velocities are ready now, save maxima
+ if(haveBins) {velocityBins->binVelSqUse(id,VelocityBins::getBodyVelSq(rb));}
+ #ifdef YADE_OPENMP
+ Real& thrMaxVSq=threadMaxVelocitySq[omp_get_thread_num()]; thrMaxVSq=max(thrMaxVSq,rb->velocity.SquaredLength());
+ #else
+ maxVelocitySq=max(maxVelocitySq,rb->velocity.SquaredLength());
+ #endif
- Vector3r axis = rb->angularVelocity;
- Real angle = axis.Normalize();
- Quaternionr q;
- q.FromAxisAngle ( axis,angle*dt );
- rb->se3.orientation = q*rb->se3.orientation;
- if(ncb->bex.getMoveRotUsed() && ncb->bex.getRot(id)!=Vector3r::ZERO){ Vector3r r(ncb->bex.getRot(id)); Real norm=r.Normalize(); q.FromAxisAngle(r,norm); rb->se3.orientation=q*rb->se3.orientation; }
- rb->se3.orientation.Normalize();
+ Vector3r axis = rb->angularVelocity;
+ Real angle = axis.Normalize();
+ Quaternionr q;
+ q.FromAxisAngle ( axis,angle*dt );
+ rb->se3.orientation = q*rb->se3.orientation;
+ if(ncb->bex.getMoveRotUsed() && ncb->bex.getRot(id)!=Vector3r::ZERO){ Vector3r r(ncb->bex.getRot(id)); Real norm=r.Normalize(); q.FromAxisAngle(r,norm); rb->se3.orientation=q*rb->se3.orientation; }
+ rb->se3.orientation.Normalize();
- rb->se3.position += rb->velocity*dt + ncb->bex.getMove(id);
+ rb->se3.position += rb->velocity*dt + ncb->bex.getMove(id);
- if(b->isClump()) static_cast<Clump*>(b.get())->moveMembers();
+ if(b->isClump()) static_cast<Clump*>(b.get())->moveMembers();
}
-
+ #ifdef YADE_OPENMP
+ FOREACH(const Real& thrMaxVSq, threadMaxVelocitySq) { maxVelocitySq=max(maxVelocitySq,thrMaxVSq); }
+ #endif
if(haveBins) velocityBins->binVelSqFinalize();
}
Modified: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp 2009-07-13 21:39:27 UTC (rev 1863)
@@ -9,6 +9,9 @@
#include<yade/core/StandAloneEngine.hpp>
#include<Wm3Vector3.h>
+#ifdef YADE_OPENMP
+ #include<omp.h>
+#endif
/*! An engine that can replace the usual series of engines used for integrating the laws of motion.
@@ -43,10 +46,17 @@
Real damping;
/// store square of max. velocity, for informative purposes; computed again at every step
Real maxVelocitySq;
+ #ifdef YADE_OPENMP
+ vector<Real> threadMaxVelocitySq;
+ #endif
/// velocity bins (not used if not created)
shared_ptr<VelocityBins> velocityBins;
virtual void action(MetaBody *);
- NewtonsDampedLaw(): damping(0.2), maxVelocitySq(-1){}
+ NewtonsDampedLaw(): damping(0.2), maxVelocitySq(-1){
+ #ifdef YADE_OPENMP
+ threadMaxVelocitySq.resize(omp_get_max_threads());
+ #endif
+ }
REGISTER_ATTRIBUTES(StandAloneEngine,(damping)(maxVelocitySq));
REGISTER_CLASS_AND_BASE(NewtonsDampedLaw,StandAloneEngine);
Modified: trunk/scripts/simple-scene.py
===================================================================
--- trunk/scripts/simple-scene.py 2009-07-13 19:31:44 UTC (rev 1862)
+++ trunk/scripts/simple-scene.py 2009-07-13 21:39:27 UTC (rev 1863)
@@ -55,21 +55,21 @@
## This MetaEngine acts on all PhysicalActions and selects the right EngineUnit base on type of the PhysicalAction.
#
# note that following 4 engines (till the end) can be replaced by an optimized monolithic version:
- # NewtonsDampedLaw(damping=0.1)
+ NewtonsDampedLaw(damping=0.1)
#
- PhysicalActionDamper([
- CundallNonViscousForceDamping(damping=0.2),
- CundallNonViscousMomentumDamping(damping=0.2)
- ]),
- ## Now we have forces and momenta acting on bodies. Newton's law calculates acceleration that corresponds to them.
- PhysicalActionApplier([
- NewtonsForceLaw(),
- NewtonsMomentumLaw(),
- ]),
- ## Acceleration results in velocity change. Integrating the velocity over dt, position of the body will change.
- PhysicalParametersMetaEngine([LeapFrogPositionIntegrator()]),
- ## Angular acceleration changes angular velocity, resulting in position and/or orientation change of the body.
- PhysicalParametersMetaEngine([LeapFrogOrientationIntegrator()])
+# PhysicalActionDamper([
+# CundallNonViscousForceDamping(damping=0.2),
+# CundallNonViscousMomentumDamping(damping=0.2)
+# ]),
+# ## Now we have forces and momenta acting on bodies. Newton's law calculates acceleration that corresponds to them.
+# PhysicalActionApplier([
+# NewtonsForceLaw(),
+# NewtonsMomentumLaw(),
+# ]),
+# ## Acceleration results in velocity change. Integrating the velocity over dt, position of the body will change.
+# PhysicalParametersMetaEngine([LeapFrogPositionIntegrator()]),
+# ## Angular acceleration changes angular velocity, resulting in position and/or orientation change of the body.
+# PhysicalParametersMetaEngine([LeapFrogOrientationIntegrator()])
]
