yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01478
[svn] r1865 - in trunk: extra pkg/dem pkg/dem/Engine/DeusExMachina pkg/dem/Engine/StandAloneEngine
Author: eudoxos
Date: 2009-07-14 23:03:51 +0200 (Tue, 14 Jul 2009)
New Revision: 1865
Added:
trunk/pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp
Removed:
trunk/extra/spherical-dem-simulator/
trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp
Modified:
trunk/extra/SConscript
trunk/pkg/dem/ConcretePM.cpp
trunk/pkg/dem/SConscript
Log:
1. Move NewtonsDampedLaw to StandAloneEngine directory
2. Add omp critical to ConcretePM
3. Remove SphericalDEMSimulator
Modified: trunk/extra/SConscript
===================================================================
--- trunk/extra/SConscript 2009-07-14 15:41:54 UTC (rev 1864)
+++ trunk/extra/SConscript 2009-07-14 21:03:51 UTC (rev 1865)
@@ -4,23 +4,6 @@
import os.path, os
env.Install('$PREFIX/lib/yade$SUFFIX/extra',[
- env.SharedLibrary('SphericalDEMSimulator',
- ['spherical-dem-simulator/PersistentAloneSAPCollider.cpp',
- 'spherical-dem-simulator/SphericalDEMSimulator.cpp'],
- LIBS=env['LIBS']+['yade-base',
- 'MacroMicroElasticRelationships',
- 'BodyMacroParameters',
- 'CundallNonViscousDamping',
- 'CundallNonViscousDamping',
- 'PhysicalActionDamper',
- 'InteractionPhysicsMetaEngine',
- 'GravityEngines',
- 'Sphere',
- 'yade-multimethods',
- 'yade-factory',
- 'yade-serialization',
- 'yade-base',
- ]),
env.SharedLibrary('Tetra',['tetra/Tetra.cpp'],LIBS=env['LIBS']+['Tetrahedron',
'ParticleParameters',
Modified: trunk/pkg/dem/ConcretePM.cpp
===================================================================
--- trunk/pkg/dem/ConcretePM.cpp 2009-07-14 15:41:54 UTC (rev 1864)
+++ trunk/pkg/dem/ConcretePM.cpp 2009-07-14 21:03:51 UTC (rev 1865)
@@ -196,7 +196,9 @@
if(isCohesive){
const shared_ptr<Body>& body1=Body::byId(I->getId1(),rootBody), body2=Body::byId(I->getId2(),rootBody); assert(body1); assert(body2);
const shared_ptr<CpmMat>& rbp1=YADE_PTR_CAST<CpmMat>(body1->physicalParameters), rbp2=YADE_PTR_CAST<CpmMat>(body2->physicalParameters);
- rbp1->numBrokenCohesive+=1; rbp2->numBrokenCohesive+=1; rbp1->epsPlBroken+=epsPlSum; rbp2->epsPlBroken+=epsPlSum;
+ // nice article about openMP::critical vs. scoped locks: http://www.thinkingparallel.com/2006/08/21/scoped-locking-vs-critical-in-openmp-a-personal-shootout/
+ #pragma omp critical
+ { rbp1->numBrokenCohesive+=1; rbp2->numBrokenCohesive+=1; rbp1->epsPlBroken+=epsPlSum; rbp2->epsPlBroken+=epsPlSum; }
}
rootBody->interactions->requestErase(I->getId1(),I->getId2());
return;
Deleted: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-07-14 15:41:54 UTC (rev 1864)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-07-14 21:03:51 UTC (rev 1865)
@@ -1,135 +0,0 @@
-/*************************************************************************
- Copyright (C) 2008 by Bruno Chareyre *
-* bruno.chareyre@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include"NewtonsDampedLaw.hpp"
-#include<yade/core/MetaBody.hpp>
-#include<yade/pkg-common/RigidBodyParameters.hpp>
-#include<yade/pkg-common/VelocityBins.hpp>
-#include<yade/lib-base/yadeWm3Extra.hpp>
-#include<yade/pkg-dem/Clump.hpp>
-
-YADE_PLUGIN("NewtonsDampedLaw");
-
-void NewtonsDampedLaw::cundallDamp(const Real& dt, const Vector3r& f, const Vector3r& velocity, Vector3r& acceleration, const Vector3r& m, const Vector3r& angularVelocity, Vector3r& angularAcceleration){
- for(int i=0; i<3; i++){
- angularAcceleration[i]*= 1 - damping*Mathr::Sign ( m[i]*(angularVelocity[i] + (Real) 0.5 *dt*angularAcceleration[i]) );
- acceleration [i]*= 1 - damping*Mathr::Sign ( f[i]*(velocity [i] + (Real) 0.5 *dt*acceleration [i]) );
- }
-}
-
-void NewtonsDampedLaw::handleClumpMember(MetaBody* ncb, const body_id_t memberId, RigidBodyParameters* clumpRBP){
- const shared_ptr<Body>& b=Body::byId(memberId,ncb);
- assert(b->isClumpMember());
- RigidBodyParameters* rb=YADE_CAST<RigidBodyParameters*>(b->physicalParameters.get());
- const Vector3r& m=ncb->bex.getTorque(memberId); const Vector3r& f=ncb->bex.getForce(memberId);
- Vector3r diffClumpAccel=f/clumpRBP->mass;
- // angular acceleration from: normal torque + torque generated by the force WRT particle centroid on the clump centroid
- Vector3r diffClumpAngularAccel=diagDiv(m,clumpRBP->inertia)+diagDiv((rb->se3.position-clumpRBP->se3.position).Cross(f),clumpRBP->inertia);
- // damp increment of accels on the clump, using velocities of the clump MEMBER
- cundallDamp(ncb->dt,f,rb->velocity,diffClumpAccel,m,rb->angularVelocity,diffClumpAngularAccel);
- // clumpRBP->{acceleration,angularAcceleration} are reset byt Clump::moveMembers, it is ok to just increment here
- clumpRBP->acceleration+=diffClumpAccel;
- clumpRBP->angularAcceleration+=diffClumpAngularAccel;
- if(haveBins) velocityBins->binVelSqUse(memberId,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
-}
-
-void NewtonsDampedLaw::action(MetaBody * ncb)
-{
- ncb->bex.sync();
- Real dt=Omega::instance().getTimeStep();
- maxVelocitySq=-1;
- haveBins=(bool)velocityBins;
- if(haveBins) velocityBins->binVelSqInitialize();
-
- #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
- // 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;
- }
-
- // 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();
-
- rb->se3.position += rb->velocity*dt + ncb->bex.getMove(id);
-
- 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();
-}
-
Deleted: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp 2009-07-14 15:41:54 UTC (rev 1864)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp 2009-07-14 21:03:51 UTC (rev 1865)
@@ -1,65 +0,0 @@
-/*************************************************************************
- Copyright (C) 2008 by Bruno Chareyre *
-* bruno.chareyre@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-#pragma once
-
-#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.
-
-This engine is faster because it uses less loops and less dispatching
-
-The result is almost the same as with :
--NewtonsForceLaw
--NewtonsMomentumLaw
--LeapFrogPositionIntegrator
--LeapFrogOrientationIntegrator
--CundallNonViscousForceDamping
--CundallNonViscousMomentumDamping
-
-...but the implementation of damping is slightly different compared to CundallNonViscousForceDamping+CundallNonViscousMomentumDamping. Here, damping is dependent on predicted (undamped) velocity at t+dt/2, while the other engines use velocity at time t.
-
-Requirements :
--All dynamic bodies must have physical parameters of type (or inheriting from) BodyMacroParameters
--Physical actions must include forces and moments
-
-NOTE: Cundall damping affected dynamic simulation! See examples/dynamic_simulation_tests
-
- */
-class RigidBodyParameters;
-class VelocityBins;
-
-class NewtonsDampedLaw : public StandAloneEngine{
- inline void cundallDamp(const Real& dt, const Vector3r& f, const Vector3r& velocity, Vector3r& acceleration, const Vector3r& m, const Vector3r& angularVelocity, Vector3r& angularAcceleration);
- void handleClumpMember(MetaBody* ncb, const body_id_t memberId, RigidBodyParameters* clumpRBP);
- bool haveBins;
- public:
- ///damping coefficient for Cundall's non viscous damping
- 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){
- #ifdef YADE_OPENMP
- threadMaxVelocitySq.resize(omp_get_max_threads());
- #endif
- }
-
- REGISTER_ATTRIBUTES(StandAloneEngine,(damping)(maxVelocitySq));
- REGISTER_CLASS_AND_BASE(NewtonsDampedLaw,StandAloneEngine);
-};
-REGISTER_SERIALIZABLE(NewtonsDampedLaw);
-
Copied: trunk/pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp (from rev 1863, trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp)
Property changes on: trunk/pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: trunk/pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp (from rev 1863, trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp)
Property changes on: trunk/pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp
___________________________________________________________________
Name: svn:mergeinfo
+
Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript 2009-07-14 15:41:54 UTC (rev 1864)
+++ trunk/pkg/dem/SConscript 2009-07-14 21:03:51 UTC (rev 1865)
@@ -950,7 +950,7 @@
env.SharedLibrary('NewtonsDampedLaw',
- ['Engine/DeusExMachina/NewtonsDampedLaw.cpp'],
+ ['Engine/StandAloneEngine/NewtonsDampedLaw.cpp'],
LIBS=env['LIBS']+['yade-serialization',
'yade-base',
'yade-multimethods',
