yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06626
[Branch ~yade-dev/yade/trunk] Rev 2626: 1. Clean up NewtonIntegrator following the discussion on the list; adds 1st order damping for asp...
------------------------------------------------------------
revno: 2626
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Wed 2010-12-22 11:55:23 +0100
message:
1. Clean up NewtonIntegrator following the discussion on the list; adds 1st order damping for aspherical particles; comments welcome.
2. Add Body.intrs() in python
modified:
core/Body.cpp
core/Body.hpp
pkg/dem/NewtonIntegrator.cpp
pkg/dem/NewtonIntegrator.hpp
--
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 'core/Body.cpp'
--- core/Body.cpp 2010-08-15 16:30:00 +0000
+++ core/Body.cpp 2010-12-22 10:55:23 +0000
@@ -1,18 +1,9 @@
-/*************************************************************************
-* Copyright (C) 2004 by Olivier Galizzi *
-* olivier.galizzi@xxxxxxx *
-* Copyright (C) 2004 by Janek Kozicki *
-* cosurgi@xxxxxxxxxx *
-* *
-* 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 "Body.hpp"
-
-#include <limits>
-#include "Scene.hpp"
-#include "Omega.hpp"
+
+#include<yade/core/Body.hpp>
+#include<limits>
+#include<yade/core/Scene.hpp>
+#include<yade/core/Omega.hpp>
+#include<yade/core/InteractionContainer.hpp>
//! This could be -1 if id_t is re-typedef'ed as `int'
const Body::id_t Body::ID_NONE=Body::id_t(-1);
@@ -20,4 +11,14 @@
const shared_ptr<Body>& Body::byId(Body::id_t _id, Scene* rb){return (*((rb?rb:Omega::instance().getScene().get())->bodies))[_id];}
const shared_ptr<Body>& Body::byId(Body::id_t _id, shared_ptr<Scene> rb){return (*(rb->bodies))[_id];}
+// return list of interactions of this particle
+python::list Body::py_intrs(){
+ Scene* scene=Omega::instance().getScene().get(); // potentially unsafe with multiple simulations
+ python::list ret;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ if(I->getId1()==id || I->getId2()==id) ret.append(I);
+ }
+ return ret;
+}
=== modified file 'core/Body.hpp'
--- core/Body.hpp 2010-12-10 15:07:56 +0000
+++ core/Body.hpp 2010-12-22 10:55:23 +0000
@@ -62,6 +62,8 @@
* (otherwise, GLViewer would depend on Clump and therefore Clump would have to go to core...) */
virtual void userForcedDisplacementRedrawHook(){return;}
+ python::list py_intrs();
+
Body::id_t getId() const {return id;};
int getGroupMask() {return groupMask; };
@@ -95,7 +97,8 @@
.def_readwrite("mask",&Body::groupMask,"Shorthand for :yref:`Body::groupMask`")
.add_property("isStandalone",&Body::isStandalone,"True if this body is neither clump, nor clump member; false otherwise.")
.add_property("isClumpMember",&Body::isClumpMember,"True if this body is clump member, false otherwise.")
- .add_property("isClump",&Body::isClump,"True if this body is clump itself, false otherwise.");
+ .add_property("isClump",&Body::isClump,"True if this body is clump itself, false otherwise.")
+ .def("intrs",&Body::py_intrs,"Return all interactions in which this body participates.")
);
};
REGISTER_SERIALIZABLE(Body);
=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp 2010-12-21 22:50:34 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2010-12-22 10:55:23 +0000
@@ -14,8 +14,9 @@
YADE_PLUGIN((NewtonIntegrator));
CREATE_LOGGER(NewtonIntegrator);
+
void NewtonIntegrator::cundallDamp(const Real& dt, const Vector3r& N, const Vector3r& V, Vector3r& A){
- for(int i=0; i<3; i++) A[i]*= 1 - damping*Mathr::Sign ( N[i]*(V[i] + (Real) 0.5 *dt*A[i]) );
+ for(int i=0; i<3; i++) A[i]*= 1 - damping*Mathr::Sign ( N[i]*(V[i] + 0.5*dt*A[i]) );
}
Vector3r NewtonIntegrator::computeAccel(const Vector3r& force, const Real& mass, int blockedDOFs){
@@ -31,26 +32,28 @@
return ret;
}
-void NewtonIntegrator::handleClumpMemberAccel(Scene* scene, const Body::id_t& memberId, State* memberState, State* clumpState){
+Vector3r NewtonIntegrator::clumpMemberAccel(const Body::id_t& memberId, const State* memberState, const State* clumpState){
const Vector3r& f=scene->forces.getForce(memberId);
- Vector3r diffClumpAccel=computeAccel(f,clumpState->mass,clumpState->blockedDOFs); // use blockedDOFs of the clump
- // damp increment of accels on the clump, using velocities of the clump MEMBER
- cundallDamp(scene->dt,f,memberState->vel,diffClumpAccel);
- clumpState->accel+=diffClumpAccel;
+ Vector3r dA=computeAccel(f,clumpState->mass,clumpState->blockedDOFs);
+ cundallDamp(scene->dt,f,memberState->vel,dA); // damp using velocities of the member
+ return dA;
}
-void NewtonIntegrator::handleClumpMemberAngAccel(Scene* scene, const Body::id_t& memberId, State* memberState, State* clumpState){
+Vector3r NewtonIntegrator::clumpMemberAngAccel(const Body::id_t& memberId, const State* memberState, const State* clumpState){
// angular acceleration from: normal torque + torque generated by the force WRT particle centroid on the clump centroid
const Vector3r& m=scene->forces.getTorque(memberId); const Vector3r& f=scene->forces.getForce(memberId);
- Vector3r diffClumpAngAccel=computeAngAccel(m+(memberState->pos-clumpState->pos).cross(f),clumpState->inertia,clumpState->blockedDOFs);
- // damp increment of accels on the clump, using velocities of the clump MEMBER
- cundallDamp(scene->dt,m,memberState->angVel,diffClumpAngAccel);
- clumpState->angAccel+=diffClumpAngAccel;
+ Vector3r dAA=computeAngAccel(m+(memberState->pos-clumpState->pos).cross(f),clumpState->inertia,clumpState->blockedDOFs);
+ cundallDamp(scene->dt,m,memberState->angVel,dAA);
+ return dAA;
}
-void NewtonIntegrator::handleClumpMemberTorque(Scene* scene, const Body::id_t& memberId, State* memberState, State* clumpState, Vector3r& M){
+Vector3r NewtonIntegrator::clumpMemberTorque(const Body::id_t& memberId, const State* memberState, const State* clumpState){
const Vector3r& m=scene->forces.getTorque(memberId); const Vector3r& f=scene->forces.getForce(memberId);
- M+=(memberState->pos-clumpState->pos).cross(f)+m;
+ Vector3r dM=(memberState->pos-clumpState->pos).cross(f)+m;
+ for(int i=0; i<3; i++) if(clumpState->blockedDOFs & State::axisDOF(i,true)) dM[i]=0.;
+ cundallDamp(0,m,memberState->angVel,dM); // damp dM; dt=0 makes the acceleration term vanish, since here it is torque rather than acceleration
+ return dM;
}
-void NewtonIntegrator::saveMaximaVelocity(Scene* scene, const Body::id_t& id, State* state){
+
+void NewtonIntegrator::saveMaximaVelocity(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.squaredNorm());
@@ -133,64 +136,48 @@
}
if (likely(b->isStandalone())){
- // translate equation
state->accel=computeAccel(f,state->mass,state->blockedDOFs);
cundallDamp(dt,f,fluctVel,state->accel);
- leapfrogTranslate(scene,state,id,dt);
- // rotate equation
+ leapfrogTranslate(state,id,dt);
// exactAsphericalRot is disabled or the body is spherical
if (likely(!exactAsphericalRot || !b->isAspherical() || b->state->inertia.minCoeff()<=0.)){
state->angAccel=computeAngAccel(m,state->inertia,state->blockedDOFs);
cundallDamp(dt,m,state->angVel,state->angAccel);
- leapfrogSphericalRotate(scene,state,id,dt);
+ leapfrogSphericalRotate(state,id,dt);
} else { // exactAsphericalRot enabled & aspherical body
- // no damping in this case
- leapfrogAsphericalRotate(scene,state,id,dt,m);
+ // damp the torque
+ Vector3r mm(m); cundallDamp(0,m,state->angVel,mm);
+ leapfrogAsphericalRotate(state,id,dt,mm);
}
} else if (b->isClump()){
- // reset acceleration of the clump itself; computed from accels on constituents
- state->accel=state->angAccel=Vector3r::Zero();
- // clump mass forces
- Vector3r dLinAccel=computeAccel(f,state->mass,state->blockedDOFs);
- cundallDamp(dt,f,fluctVel,dLinAccel);
- state->accel+=dLinAccel;
- Vector3r M(m); // torque on the clump itself, will accumulate from members
- // sum force on clump memebrs
- // exactAsphericalRot enabled and clump is aspherical
- if (exactAsphericalRot && b->isAspherical() && b->state->inertia.maxCoeff()>0){
- FOREACH(Clump::MemberMap::value_type mm, static_cast<Clump*>(b->shape.get())->members){
- const Body::id_t& memberId=mm.first;
- const shared_ptr<Body>& member=Body::byId(memberId,scene); assert(member->isClumpMember());
- State* memberState=member->state.get();
- handleClumpMemberAccel(scene,memberId,memberState,state);
- handleClumpMemberTorque(scene,memberId,memberState,state,M);
- }
- // motion
- leapfrogTranslate(scene,state,id,dt);
- leapfrogAsphericalRotate(scene,state,id,dt,M);
- } else { // exactAsphericalRot disabled or clump is spherical
- Vector3r dAngAccel=computeAngAccel(M,state->inertia,state->blockedDOFs);
- cundallDamp(dt,M,state->angVel,dAngAccel);
- state->angAccel+=dAngAccel;
- FOREACH(Clump::MemberMap::value_type mm, static_cast<Clump*>(b->shape.get())->members){
- const Body::id_t& memberId=mm.first;
- const shared_ptr<Body>& member=Body::byId(memberId,scene); assert(member->isClumpMember());
- State* memberState=member->state.get();
- handleClumpMemberAccel(scene,memberId,memberState,state);
- handleClumpMemberAngAccel(scene,memberId,memberState,state);
- }
- // motion
- leapfrogTranslate(scene,state,id,dt);
- leapfrogSphericalRotate(scene,state,id,dt);
+ // aspherical if allowed, if the body is aspherical, has no zero inertia term (->NaN) and has at least one rotational DOF unblocked
+ const bool aspherical=(exactAsphericalRot && b->isAspherical() && b->state->inertia.maxCoeff()>0 && ((state->blockedDOFs&State::DOF_RXRYRZ)!=State::DOF_RXRYRZ));
+ Vector3r M(m); // torque on the clump itself, aspherical only
+ state->accel=computeAccel(f,state->mass,state->blockedDOFs); cundallDamp(dt,f,fluctVel,state->accel);
+ if(!aspherical){ state->angAccel=computeAngAccel(M,state->inertia,state->blockedDOFs); cundallDamp(dt,m,state->angVel,state->angAccel); }
+ if(state->blockedDOFs!=State::DOF_ALL){ // if the clump has all DoFs blocked, forces on members would be computed uselessly
+ FOREACH(Clump::MemberMap::value_type mm, static_cast<Clump*>(b->shape.get())->members){
+ const Body::id_t& memberId=mm.first;
+ const shared_ptr<Body>& member=Body::byId(memberId,scene); assert(member->isClumpMember());
+ State* memberState=member->state.get();
+ state->accel+=clumpMemberAccel(memberId,memberState,state);
+ if(aspherical) M+=clumpMemberTorque(memberId,memberState,state);
+ else state->angAccel+=clumpMemberAngAccel(memberId,memberState,state);
+ }
}
+ // translation
+ leapfrogTranslate(state,id,dt);
+ // rotation
+ if(aspherical) leapfrogAsphericalRotate(state,id,dt,M);
+ else leapfrogSphericalRotate(state,id,dt);
Clump::moveMembers(b,scene);
// save max velocity for all members of the clump
FOREACH(Clump::MemberMap::value_type mm, static_cast<Clump*>(b->shape.get())->members){
const shared_ptr<Body>& member=Body::byId(mm.first,scene); assert(member->isClumpMember());
- saveMaximaVelocity(scene,mm.first,member->state.get());
+ saveMaximaVelocity(mm.first,member->state.get());
}
}
- saveMaximaVelocity(scene,id,state);
+ saveMaximaVelocity(id,state);
// process callbacks
for(size_t i=0; i<callbacksSize; i++){
cerr<<"<"<<b->id<<",cb="<<callbacks[i]<<",scene="<<callbacks[i]->scene<<">"; // <<",force="<<callbacks[i]->scene->forces.getForce(b->id)<<">";
@@ -204,7 +191,7 @@
if(scene->isPeriodic) { prevCellSize=scene->cell->getSize(); prevVelGrad=scene->cell->velGrad; }
}
-inline void NewtonIntegrator::leapfrogTranslate(Scene* scene, State* state, const Body::id_t& id, const Real& dt){
+void NewtonIntegrator::leapfrogTranslate(State* state, const Body::id_t& id, const Real& dt){
if (scene->forces.getMoveRotUsed()) state->pos+=scene->forces.getMove(id);
if (homoDeform==Cell::HOMO_VEL || homoDeform==Cell::HOMO_VEL_2ND) {
// update velocity reflecting changes in the macroscopic velocity field, making the problem homothetic.
@@ -223,7 +210,7 @@
state->pos+=state->vel*dt;
}
-inline void NewtonIntegrator::leapfrogSphericalRotate(Scene* scene, State* state, const Body::id_t& id, const Real& dt )
+void NewtonIntegrator::leapfrogSphericalRotate(State* state, const Body::id_t& id, const Real& dt )
{
state->angVel+=dt*state->angAccel;
Vector3r axis = state->angVel;
@@ -242,15 +229,7 @@
state->ori.normalize();
}
-void NewtonIntegrator::leapfrogAsphericalRotate(Scene* scene, State* state, const Body::id_t& id, const Real& dt, const Vector3r& M){
- // if all DoFs are blocked then return; warn if there is angular acceleration, since it will be ignored
- if((state->blockedDOFs & State::DOF_RXRYRZ)==State::DOF_RXRYRZ){
- /* necessarily angular acceleration is zero, we can call the spherical integrator safely */
- if(state->angVel!=Vector3r::Zero()) leapfrogSphericalRotate(scene,state,id,dt); /* in most cases, angVel will be zero; no need to integrate anything, just return */
- return;
- }
- // if only some rotations are allowed, warn, since it is not handled (yet?)
- if((state->blockedDOFs & State::DOF_RXRYRZ)!=0) LOG_WARN("Aspherical body #"<<id<<" with some rotational DoFs blocked is not supported (blocked DoFs ignored).");
+void NewtonIntegrator::leapfrogAsphericalRotate(State* state, const Body::id_t& id, const Real& dt, const Vector3r& M){
Matrix3r A=state->ori.conjugate().toRotationMatrix(); // 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/NewtonIntegrator.hpp'
--- pkg/dem/NewtonIntegrator.hpp 2010-12-21 15:31:40 +0000
+++ pkg/dem/NewtonIntegrator.hpp 2010-12-22 10:55:23 +0000
@@ -23,19 +23,20 @@
class NewtonIntegrator : public GlobalEngine{
inline void cundallDamp(const Real& dt, const Vector3r& N, const Vector3r& V, Vector3r& A);
- inline void handleClumpMemberAccel(Scene* ncb, const Body::id_t& memberId, State* memberState, State* clumpState);
- inline void handleClumpMemberAngAccel(Scene* ncb, const Body::id_t& memberId, State* memberState, State* clumpState);
- inline void handleClumpMemberTorque(Scene* ncb, const Body::id_t& memberId, State* memberState, State* clumpState, Vector3r& M);
- inline void saveMaximaVelocity(Scene* ncb, const Body::id_t& id, State* state);
+ inline void saveMaximaVelocity(const Body::id_t& id, State* state);
bool haveBins;
- inline void leapfrogTranslate(Scene* ncb, State* state, const Body::id_t& id, const Real& dt); // leap-frog translate
- inline void leapfrogSphericalRotate(Scene* ncb, State* state, const Body::id_t& id, const Real& dt); // leap-frog rotate of spherical body
- inline void leapfrogAsphericalRotate(Scene* ncb, State* state, const Body::id_t& id, const Real& dt, const Vector3r& M); // leap-frog rotate of aspherical body
+ inline void leapfrogTranslate(State*, const Body::id_t& id, const Real& dt); // leap-frog translate
+ inline void leapfrogSphericalRotate(State*, const Body::id_t& id, const Real& dt); // leap-frog rotate of spherical body
+ inline void leapfrogAsphericalRotate(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);
// compute linear and angular acceleration, respecting State::blockedDOFs
Vector3r computeAccel(const Vector3r& force, const Real& mass, int blockedDOFs);
Vector3r computeAngAccel(const Vector3r& torque, const Vector3r& inertia, int blockedDOFs);
+ // compute contribution of clump member
+ Vector3r clumpMemberAccel(const Body::id_t&, const State* clumpState, const State* memberState);
+ Vector3r clumpMemberAngAccel(const Body::id_t&, const State* clumpState, const State* memberState);
+ Vector3r clumpMemberTorque(const Body::id_t&, const State* clumpState, const State* memberState);
// whether the cell has changed from the previous step
bool cellChanged;
Follow ups
