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