yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02456
[Branch ~yade-dev/yade/trunk] Rev 1840: 1. NewtonDampedLaw: rename functions and remove duplicate code as Vaclav suggest (thanks!)
Merge authors:
Sergei D. <sega@think>
------------------------------------------------------------
revno: 1840 [merge]
committer: Sergei D. <sega@think>
branch nick: trunk
timestamp: Wed 2009-12-02 12:26:16 +0300
message:
1. NewtonDampedLaw: rename functions and remove duplicate code as Vaclav suggest (thanks!)
2. Add spheres ids to VTKRecorder
modified:
pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp
pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp
pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp
pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp
scripts/test/clump-hopper-viscoelastic.py
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.
=== modified file 'pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp'
--- pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp 2009-12-01 14:56:39 +0000
+++ pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp 2009-12-02 09:26:16 +0000
@@ -31,22 +31,27 @@
if((blockedDOFs & State::DOF_RY)!=0) v[1]=0;
if((blockedDOFs & State::DOF_RZ)!=0) v[2]=0;
}
-void NewtonsDampedLaw::handleClumpMember(World* ncb, const body_id_t memberId, State* clumpState){
- const shared_ptr<Body>& b=Body::byId(memberId,ncb);
- assert(b->isClumpMember());
- State* state=b->state.get();
- const Vector3r& m=ncb->bex.getTorque(memberId); const Vector3r& f=ncb->bex.getForce(memberId);
+void NewtonsDampedLaw::handleClumpMemberAccel(World* ncb, const body_id_t& memberId, State* memberState, State* clumpState){
+ const Vector3r& f=ncb->bex.getForce(memberId);
Vector3r diffClumpAccel=f/clumpState->mass;
+ // damp increment of accels on the clump, using velocities of the clump MEMBER
+ cundallDamp(ncb->dt,f,memberState->vel,diffClumpAccel);
+ clumpState->accel+=diffClumpAccel;
+}
+void NewtonsDampedLaw::handleClumpMemberAngAccel(World* ncb, const body_id_t& memberId, State* memberState, State* clumpState){
// angular acceleration from: normal torque + torque generated by the force WRT particle centroid on the clump centroid
- Vector3r diffClumpAngularAccel=diagDiv(m,clumpState->inertia)+diagDiv((state->pos-clumpState->pos).Cross(f),clumpState->inertia);
+ const Vector3r& m=ncb->bex.getTorque(memberId); const Vector3r& f=ncb->bex.getForce(memberId);
+ Vector3r diffClumpAngularAccel=diagDiv(m,clumpState->inertia)+diagDiv((memberState->pos-clumpState->pos).Cross(f),clumpState->inertia);
// damp increment of accels on the clump, using velocities of the clump MEMBER
- //cundallDamp(ncb->dt,f,state->vel,diffClumpAccel,m,state->angVel,diffClumpAngularAccel);
- cundallDamp(ncb->dt,f,state->vel,diffClumpAccel);
- cundallDamp(ncb->dt,m,state->angVel,diffClumpAngularAccel);
- // clumpState->{acceleration,angularAcceleration} are reset byt Clump::moveMembers, it is ok to just increment here
- clumpState->accel+=diffClumpAccel;
+ cundallDamp(ncb->dt,m,memberState->angVel,diffClumpAngularAccel);
clumpState->angAccel+=diffClumpAngularAccel;
- if(haveBins) velocityBins->binVelSqUse(memberId,VelocityBins::getBodyVelSq(state));
+}
+void NewtonsDampedLaw::handleClumpMemberTorque(World* ncb, const body_id_t& memberId, State* memberState, State* clumpState, Vector3r& M){
+ const Vector3r& m=ncb->bex.getTorque(memberId); const Vector3r& f=ncb->bex.getForce(memberId);
+ M+=(memberState->pos-clumpState->pos).Cross(f)+m;
+}
+void NewtonsDampedLaw::saveMaximaVelocity(World* ncb, const body_id_t& id, State* state){
+ if(haveBins) velocityBins->binVelSqUse(id,VelocityBins::getBodyVelSq(state));
#ifdef YADE_OPENMP
Real& thrMaxVSq=threadMaxVelocitySq[omp_get_thread_num()]; thrMaxVSq=max(thrMaxVSq,state->vel.SquaredLength());
#else
@@ -77,13 +82,7 @@
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(state));}
- #ifdef YADE_OPENMP
- Real& thrMaxVSq=threadMaxVelocitySq[omp_get_thread_num()]; thrMaxVSq=max(thrMaxVSq,state->vel.SquaredLength());
- #else
- maxVelocitySq=max(maxVelocitySq,state->vel.SquaredLength());
- #endif
+ saveMaximaVelocity(ncb,id,state);
continue;
}
@@ -92,76 +91,57 @@
const Vector3r& f=ncb->bex.getForce(id);
state->accel=f/state->mass;
cundallDamp(dt,f,state->vel,state->accel);
- lfTranslate(ncb,state,id,dt);
+ leapfrogTranslate(ncb,state,id,dt);
// rotate equation
- if (/*b->isSpheral || */ !accRigidBodyRot){ // spheral body or accRigidBodyRot disabled
+ if (/*b->isSpheral || */ !exactAsphericalRot){ // spheral body or exactAsphericalRot disabled
const Vector3r& m=ncb->bex.getTorque(id);
state->angAccel=diagDiv(m,state->inertia);
cundallDamp(dt,m,state->angVel,state->angAccel);
- lfSpheralRotate(ncb,state,id,dt);
- } else { // non spheral body and accRigidBodyRot enabled
+ leapfrogSphericalRotate(ncb,state,id,dt);
+ } else { // non spheral body and exactAsphericalRot enabled
const Vector3r& m=ncb->bex.getTorque(id);
- lfRigidBodyRotate(ncb,state,id,dt,m);
+ leapfrogAsphericalRotate(ncb,state,id,dt,m);
}
} else if (b->isClump()){
- state->accel=state->angAccel=Vector3r::ZERO; // to make sure; should be reset in Clump::moveMembers
- if (accRigidBodyRot){
- // forces applied to clump proper, if there are such
- const Vector3r& f=ncb->bex.getForce(id);
- Vector3r dLinAccel=f/state->mass;
- cundallDamp(dt,f,state->vel,dLinAccel);
- state->accel+=dLinAccel;
- const Vector3r& m=ncb->bex.getTorque(id);
- Vector3r M(m);
- // sum forces on clump members
- FOREACH(Clump::memberMap::value_type mm, static_cast<Clump*>(b.get())->members){
- const body_id_t memberId=mm.first;
- const shared_ptr<Body>& member=Body::byId(memberId,ncb);
- assert(member->isClumpMember());
- State* memberState=member->state.get();
- // Linear acceleration
- const Vector3r& f=ncb->bex.getForce(memberId);
- Vector3r diffClumpAccel=f/state->mass;
- // damp increment of accel on the clump, using velocities of the clump MEMBER
- cundallDamp(dt,f,memberState->vel,diffClumpAccel);
- state->accel+=diffClumpAccel;
- // Momentum
- const Vector3r& m=ncb->bex.getTorque(memberId);
- M+=(memberState->pos-state->pos).Cross(f)+m;
- if(haveBins) velocityBins->binVelSqUse(memberId,VelocityBins::getBodyVelSq(memberState));
- #ifdef YADE_OPENMP
- Real& thrMaxVSq=threadMaxVelocitySq[omp_get_thread_num()]; thrMaxVSq=max(thrMaxVSq,memberState->vel.SquaredLength());
- #else
- maxVelocitySq=max(maxVelocitySq,memberState->vel.SquaredLength());
- #endif
- }
- // motion
- lfTranslate(ncb,state,id,dt);
- lfRigidBodyRotate(ncb,state,id,dt,M);
- } else { // accRigidBodyRot disabled
- // 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,state);
- }
- // forces applied to clump proper, if there are such
- const Vector3r& m=ncb->bex.getTorque(id); const Vector3r& f=ncb->bex.getForce(id);
- Vector3r dLinAccel=f/state->mass, dAngAccel=diagDiv(m,state->inertia);
- cundallDamp(dt,f,state->vel,dLinAccel); cundallDamp(dt,m,state->angVel,dAngAccel);
- state->accel+=dLinAccel; state->angAccel+=dAngAccel;
- // motion
- lfTranslate(ncb,state,id,dt);
- lfSpheralRotate(ncb,state,id,dt);
+ // clump mass forces
+ const Vector3r& f=ncb->bex.getForce(id);
+ Vector3r dLinAccel=f/state->mass;
+ cundallDamp(dt,f,state->vel,dLinAccel);
+ state->accel+=dLinAccel;
+ const Vector3r& m=ncb->bex.getTorque(id);
+ Vector3r M(m);
+ // sum force on clump memebrs
+ if (exactAsphericalRot){
+ FOREACH(Clump::memberMap::value_type mm, static_cast<Clump*>(b.get())->members){
+ const body_id_t& memberId=mm.first;
+ const shared_ptr<Body>& member=Body::byId(memberId,ncb); assert(b->isClumpMember());
+ State* memberState=member->state.get();
+ handleClumpMemberAccel(ncb,memberId,memberState,state);
+ handleClumpMemberTorque(ncb,memberId,memberState,state,M);
+ saveMaximaVelocity(ncb,memberId,memberState);
+ }
+ // motion
+ leapfrogTranslate(ncb,state,id,dt);
+ leapfrogAsphericalRotate(ncb,state,id,dt,M);
+ } else { // exactAsphericalRot disabled
+ Vector3r dAngAccel=diagDiv(M,state->inertia);
+ cundallDamp(dt,M,state->angVel,dAngAccel);
+ state->angAccel+=dAngAccel;
+ FOREACH(Clump::memberMap::value_type mm, static_cast<Clump*>(b.get())->members){
+ const body_id_t& memberId=mm.first;
+ const shared_ptr<Body>& member=Body::byId(memberId,ncb); assert(b->isClumpMember());
+ State* memberState=member->state.get();
+ handleClumpMemberAccel(ncb,memberId,memberState,state);
+ handleClumpMemberAngAccel(ncb,memberId,memberState,state);
+ saveMaximaVelocity(ncb,memberId,memberState);
+ }
+ // motion
+ leapfrogTranslate(ncb,state,id,dt);
+ leapfrogSphericalRotate(ncb,state,id,dt);
}
static_cast<Clump*>(b.get())->moveMembers();
}
-
- // save maxima velocity
- if(haveBins) {velocityBins->binVelSqUse(id,VelocityBins::getBodyVelSq(state));}
- #ifdef YADE_OPENMP
- Real& thrMaxVSq=threadMaxVelocitySq[omp_get_thread_num()]; thrMaxVSq=max(thrMaxVSq,state->vel.SquaredLength());
- #else
- maxVelocitySq=max(maxVelocitySq,state->vel.SquaredLength());
- #endif
+ saveMaximaVelocity(ncb,id,state);
}
#ifdef YADE_OPENMP
FOREACH(const Real& thrMaxVSq, threadMaxVelocitySq) { maxVelocitySq=max(maxVelocitySq,thrMaxVSq); }
@@ -169,13 +149,13 @@
if(haveBins) velocityBins->binVelSqFinalize();
}
-inline void NewtonsDampedLaw::lfTranslate(World* ncb, State* state, const body_id_t& id, const Real& dt )
+inline void NewtonsDampedLaw::leapfrogTranslate(World* ncb, State* state, const body_id_t& id, const Real& dt )
{
blockTranslateDOFs(state->blockedDOFs, state->accel);
state->vel+=dt*state->accel;
state->pos += state->vel*dt + ncb->bex.getMove(id);
}
-inline void NewtonsDampedLaw::lfSpheralRotate(World* ncb, State* state, const body_id_t& id, const Real& dt )
+inline void NewtonsDampedLaw::leapfrogSphericalRotate(World* ncb, State* state, const body_id_t& id, const Real& dt )
{
blockRotateDOFs(state->blockedDOFs, state->angAccel);
state->angVel+=dt*state->angAccel;
@@ -185,7 +165,7 @@
if(ncb->bex.getMoveRotUsed() && ncb->bex.getRot(id)!=Vector3r::ZERO){ Vector3r r(ncb->bex.getRot(id)); Real norm=r.Normalize(); Quaternionr q; q.FromAxisAngle(r,norm); state->ori=q*state->ori; }
state->ori.Normalize();
}
-void NewtonsDampedLaw::lfRigidBodyRotate(World* ncb, State* state, const body_id_t& id, const Real& dt, const Vector3r& M){
+void NewtonsDampedLaw::leapfrogAsphericalRotate(World* ncb, State* state, const body_id_t& id, const Real& dt, const Vector3r& M){
Matrix3r A; state->ori.Conjugate().ToRotationMatrix(A); // rotation matrix from global to local r.f.
const Vector3r l_n = state->angMom + dt/2 * M; // global angular momentum at time n
const Vector3r l_b_n = A*l_n; // local angular momentum at time n
=== modified file 'pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp'
--- pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp 2009-12-01 14:56:39 +0000
+++ pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp 2009-12-02 09:26:16 +0000
@@ -39,11 +39,14 @@
class NewtonsDampedLaw : public StandAloneEngine{
inline void cundallDamp(const Real& dt, const Vector3r& N, const Vector3r& V, Vector3r& A);
- void handleClumpMember(World* ncb, const body_id_t memberId, State* clumpState);
+ inline void handleClumpMemberAccel(World* ncb, const body_id_t& memberId, State* memberState, State* clumpState);
+ inline void handleClumpMemberAngAccel(World* ncb, const body_id_t& memberId, State* memberState, State* clumpState);
+ inline void handleClumpMemberTorque(World* ncb, const body_id_t& memberId, State* memberState, State* clumpState, Vector3r& M);
+ inline void saveMaximaVelocity(World* ncb, const body_id_t& id, State* state);
bool haveBins;
- inline void lfTranslate(World* ncb, State* state, const body_id_t& id, const Real& dt); // leap-frog translate
- inline void lfSpheralRotate(World* ncb, State* state, const body_id_t& id, const Real& dt); // leap-frog rotate of spheral body
- inline void lfRigidBodyRotate(World* ncb, State* state, const body_id_t& id, const Real& dt, const Vector3r& M); // leap-frog rotate of rigid (non symmetric) body
+ inline void leapfrogTranslate(World* ncb, State* state, const body_id_t& id, const Real& dt); // leap-frog translate
+ inline void leapfrogSphericalRotate(World* ncb, State* state, const body_id_t& id, const Real& dt); // leap-frog rotate of spherical body
+ inline void leapfrogAsphericalRotate(World* ncb, State* state, const body_id_t& id, const Real& dt, const Vector3r& M); // leap-frog rotate of aspherical body
Quaternionr DotQ(const Vector3r& angVel, const Quaternionr& Q);
inline void blockTranslateDOFs(unsigned blockedDOFs, Vector3r& v);
inline void blockRotateDOFs(unsigned blockedDOFs, Vector3r& v);
@@ -53,21 +56,21 @@
Real damping;
/// store square of max. velocity, for informative purposes; computed again at every step
Real maxVelocitySq;
- /// Enable of the accurate rigid body rotation integrator
- bool accRigidBodyRot;
+ /// Enable of the exact aspherical body rotation integrator
+ bool exactAsphericalRot;
#ifdef YADE_OPENMP
vector<Real> threadMaxVelocitySq;
#endif
/// velocity bins (not used if not created)
shared_ptr<VelocityBins> velocityBins;
virtual void action(World *);
- NewtonsDampedLaw(): damping(0.2), maxVelocitySq(-1), accRigidBodyRot(false){
+ NewtonsDampedLaw(): damping(0.2), maxVelocitySq(-1), exactAsphericalRot(false){
#ifdef YADE_OPENMP
threadMaxVelocitySq.resize(omp_get_max_threads());
#endif
}
- REGISTER_ATTRIBUTES(StandAloneEngine,(damping)(maxVelocitySq)(accRigidBodyRot));
+ REGISTER_ATTRIBUTES(StandAloneEngine,(damping)(maxVelocitySq)(exactAsphericalRot));
REGISTER_CLASS_AND_BASE(NewtonsDampedLaw,StandAloneEngine);
DECLARE_LOGGER;
};
=== modified file 'pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp'
--- pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-12-01 14:56:39 +0000
+++ pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-12-02 09:26:16 +0000
@@ -50,7 +50,8 @@
else if(rec=="colors") recActive[REC_COLORS]=true;
else if(rec=="cpm") recActive[REC_CPM]=true;
else if(rec=="intr") recActive[REC_INTR]=true;
- else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: spheres, velocity, facets, colors, cpm, intr). Ignored.");
+ else if(rec=="ids") recActive[REC_IDS]=true;
+ else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: spheres, velocity, facets, colors, cpm, intr, ids). Ignored.");
}
// cpm needs interactions
if(recActive[REC_CPM]) recActive[REC_INTR]=true;
@@ -61,6 +62,9 @@
vtkSmartPointer<vtkFloatArray> radii = vtkSmartPointer<vtkFloatArray>::New();
radii->SetNumberOfComponents(1);
radii->SetName("radii");
+ vtkSmartPointer<vtkFloatArray> spheresIds = vtkSmartPointer<vtkFloatArray>::New();
+ spheresIds->SetNumberOfComponents(1);
+ spheresIds->SetName("IDS");
vtkSmartPointer<vtkFloatArray> spheresColors = vtkSmartPointer<vtkFloatArray>::New();
spheresColors->SetNumberOfComponents(3);
spheresColors->SetName("colors");
@@ -143,6 +147,7 @@
pid[0] = spheresPos->InsertNextPoint(pos[0], pos[1], pos[2]);
spheresCells->InsertNextCell(1,pid);
radii->InsertNextValue(sphere->radius);
+ if (recActive[REC_IDS]) spheresIds->InsertNextValue(b->getId());
if (recActive[REC_COLORS])
{
const Vector3r& color = sphere->diffuseColor;
@@ -210,6 +215,7 @@
spheresUg->SetPoints(spheresPos);
spheresUg->SetCells(VTK_VERTEX, spheresCells);
spheresUg->GetPointData()->AddArray(radii);
+ if (recActive[REC_IDS]) spheresUg->GetPointData()->AddArray(spheresIds);
if (recActive[REC_COLORS]) spheresUg->GetPointData()->AddArray(spheresColors);
if (recActive[REC_VELOCITY]) {
spheresUg->GetPointData()->AddArray(spheresVelocity);
=== modified file 'pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp'
--- pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp 2009-12-01 14:56:39 +0000
+++ pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp 2009-12-02 09:26:16 +0000
@@ -4,7 +4,7 @@
class VTKRecorder: public PeriodicEngine {
public:
- enum {REC_SPHERES=0,REC_FACETS,REC_COLORS,REC_CPM,REC_INTR,REC_VELOCITY,REC_SENTINEL};
+ enum {REC_SPHERES=0,REC_FACETS,REC_COLORS,REC_CPM,REC_INTR,REC_VELOCITY,REC_IDS,REC_SENTINEL};
//! A stuff to record: spheres,facets,colors
vector<string> recorders;
string fileName;
=== modified file 'scripts/test/clump-hopper-viscoelastic.py'
--- scripts/test/clump-hopper-viscoelastic.py 2009-11-30 12:50:41 +0000
+++ scripts/test/clump-hopper-viscoelastic.py 2009-12-02 09:26:16 +0000
@@ -47,8 +47,9 @@
tblIds=O.bodies.append(pack.gtsSurface2Facets(table.faces(),material=facetMat,color=(0,1,0)))
# Create clumps...
+clumpColor=(0.0, 0.5, 0.5)
for k,l in itertools.product(arange(0,10),arange(0,10)):
- clpId,sphId=O.bodies.appendClumped([utils.sphere(Vector3(x0t+Rs*(k*4+2),y0t+Rs*(l*4+2),i*Rs*2+zt),Rs,color=(0,0,1),material=dfltSpheresMat) for i in xrange(4)])
+ clpId,sphId=O.bodies.appendClumped([utils.sphere(Vector3(x0t+Rs*(k*4+2),y0t+Rs*(l*4+2),i*Rs*2+zt),Rs,color=clumpColor,material=dfltSpheresMat) for i in xrange(4)])
for id in sphId:
s=O.bodies[id]
p=utils.getViscoelasticFromSpheresInteraction(s.state['mass'],tc,en,es)
@@ -57,8 +58,9 @@
#O.bodies[clpId].state.blockedDOFs=['x','y']
# ... and spheres
+spheresColor=(0.4, 0.4, 0.4)
for k,l in itertools.product(arange(0,9),arange(0,9)):
- sphAloneId=O.bodies.append( [utils.sphere( Vector3(x0t+Rs*(k*4+4),y0t+Rs*(l*4+4),i*Rs*2.3+zt),Rs,color=(0,1,0),material=dfltSpheresMat) for i in xrange(4) ] )
+ sphAloneId=O.bodies.append( [utils.sphere( Vector3(x0t+Rs*(k*4+4),y0t+Rs*(l*4+4),i*Rs*2.3+zt),Rs,color=spheresColor,material=dfltSpheresMat) for i in xrange(4) ] )
for id in sphAloneId:
s=O.bodies[id]
p=utils.getViscoelasticFromSpheresInteraction(s.state['mass'],tc,en,es)
@@ -70,14 +72,14 @@
O.engines=[
BexResetter(),
BoundingVolumeMetaEngine([InteractingSphere2AABB(),InteractingFacet2AABB()]),
- InsertionSortCollider(),
+ InsertionSortCollider(nBins=5,sweepLength=.1*Rs),
InteractionDispatchers(
[InteractingSphere2InteractingSphere4SpheresContactGeometry(), InteractingFacet2InteractingSphere4SpheresContactGeometry()],
[Ip2_SimleViscoelasticMat_SimpleViscoelasticMat_SimpleViscoelasticPhys()],
[Law2_Spheres_Viscoelastic_SimpleViscoelastic()],
),
GravityEngine(gravity=[0,0,-9.81]),
- NewtonsDampedLaw(damping=0,accRigidBodyRot=True),
+ NewtonsDampedLaw(damping=0,exactAsphericalRot=True),
#VTKRecorder(virtPeriod=0.01,fileName='/tmp/',recorders=['spheres','velocity','facets'])
]
