yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02339
[Branch ~yade-dev/yade/trunk] Rev 1812: Add to NewtonDampedLaw a modified leap-frog algorithm for accurate integration of rotation equati...
------------------------------------------------------------
revno: 1812
committer: Sergei D. <sega@think>
branch nick: clumps
timestamp: Thu 2009-11-26 13:34:14 +0300
message:
Add to NewtonDampedLaw a modified leap-frog algorithm for accurate integration of rotation equation. This is intermideate commit: new algo seems to work, but needs some coupling with previous code.
modified:
core/State.hpp
pkg/dem/DataClass/Clump.cpp
pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp
pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp
pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp
pkg/dem/Engine/StandAloneEngine/VTKRecorder.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/State.hpp'
--- core/State.hpp 2009-11-21 16:46:58 +0000
+++ core/State.hpp 2009-11-26 10:34:14 +0000
@@ -30,6 +30,8 @@
Vector3r angAccel;
Vector3r inertia;
+ Vector3r prevAngMom;
+
/// reference values
Vector3r refPos;
Quaternionr refOri;
@@ -58,9 +60,9 @@
//! Setter of blockedDOFs from list of strings (['x','rx','rz'] → DOF_X | DOR_RX | DOF_RZ)
void blockedDOFs_vec_set(const std::vector<std::string>& dofs);
- State(): se3(Vector3r::ZERO,Quaternionr::IDENTITY),pos(se3.position),vel(Vector3r::ZERO),accel(Vector3r::ZERO),mass(0.),ori(se3.orientation),angVel(Vector3r::ZERO),angAccel(Vector3r::ZERO),inertia(Vector3r::ZERO),refPos(Vector3r::ZERO),refOri(Quaternionr::IDENTITY),blockedDOFs(DOF_NONE){}
+ State(): se3(Vector3r::ZERO,Quaternionr::IDENTITY),pos(se3.position),vel(Vector3r::ZERO),accel(Vector3r::ZERO),mass(0.),ori(se3.orientation),angVel(Vector3r::ZERO),angAccel(Vector3r::ZERO),inertia(Vector3r::ZERO),refPos(Vector3r::ZERO),refOri(Quaternionr::IDENTITY),blockedDOFs(DOF_NONE),prevAngMom(Vector3r::ZERO){}
REGISTER_CLASS_AND_BASE(State,Serializable);
- REGISTER_ATTRIBUTES(Serializable,(pos)(vel)(accel)(mass)(ori)(angVel)(angAccel)(inertia)(refPos)(refOri)(blockedDOFs));
+ REGISTER_ATTRIBUTES(Serializable,(pos)(vel)(accel)(mass)(ori)(angVel)(angAccel)(inertia)(refPos)(refOri)(blockedDOFs)(prevAngMom));
};
REGISTER_SERIALIZABLE(State);
=== modified file 'pkg/dem/DataClass/Clump.cpp'
--- pkg/dem/DataClass/Clump.cpp 2009-11-21 11:21:54 +0000
+++ pkg/dem/DataClass/Clump.cpp 2009-11-26 10:34:14 +0000
@@ -115,7 +115,8 @@
//LOG_TRACE("Clump #"<<getId()<<" moved #"<<I->first<<".");
//! FIXME: we set velocity because of damping here; but since positions are integrated after all forces applied, these velocities will be used in the NEXT step for CundallNonViscousDamping. Does that matter?!
- subState->vel=state->vel+state->angVel.Cross(I->second.position);
+ //subState->vel=state->vel+state->angVel.Cross(I->second.position);
+ subState->vel=state->vel+state->angVel.Cross(subState->pos-state->pos);
subState->angVel=state->angVel;
}
/* @bug Temporarily we reset acceleration and angularAcceleration of the clump here;
=== modified file 'pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp'
--- pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp 2009-11-23 15:36:10 +0000
+++ pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp 2009-11-26 10:34:14 +0000
@@ -16,6 +16,7 @@
#include<yade/lib-base/yadeWm3Extra.hpp>
YADE_PLUGIN((NewtonsDampedLaw));
+CREATE_LOGGER(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]) );
@@ -83,6 +84,11 @@
cundallDamp(dt,f,state->vel,state->accel,m,state->angVel,state->angAccel);
}
else if (b->isClump()){
+ if (accRigidBodyRot) {
+ accurateRigidBodyRotationIntegrator(ncb,b);
+ continue;
+ }
+
state->accel=state->angAccel=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){
@@ -136,3 +142,65 @@
#endif
if(haveBins) velocityBins->binVelSqFinalize();
}
+
+void NewtonsDampedLaw::accurateRigidBodyRotationIntegrator(MetaBody* ncb, const shared_ptr<Body>& rb){
+ Real dt=Omega::instance().getTimeStep();
+ State* state=rb->state.get();
+ const body_id_t id=rb->getId();
+
+ state->accel=state->angAccel=Vector3r::ZERO; // to make sure; should be reset in Clump::moveMembers
+
+ // sum of forces and torques
+ const Vector3r& f=ncb->bex.getForce(id); const Vector3r& m=ncb->bex.getTorque(id);
+ Vector3r F(f), M(m);
+ FOREACH(Clump::memberMap::value_type mm, static_cast<Clump*>(rb.get())->members){
+ const body_id_t memberId=mm.first;
+ const shared_ptr<Body>& b=Body::byId(memberId,ncb);
+ assert(b->isClumpMember());
+ State* bs=b->state.get();
+ const Vector3r& f=ncb->bex.getForce(memberId); const Vector3r& m=ncb->bex.getTorque(memberId);
+ F+=f; M+=(bs->pos-state->pos).Cross(f)+m;
+ }
+ // translation equation
+ state->accel=F/state->mass;
+ state->vel+=dt*state->accel;
+ state->pos+=state->vel*dt;//+ncb->bex.getMove(id);
+ // rotation equation
+ Vector3r l_n = state->prevAngMom + dt/2 * M; // global angular momentum at time n
+ Vector3r l_b_n = state->ori.Conjugate().Rotate(l_n); // local angular momentum at time n
+ Vector3r angVel_b_n = diagDiv(l_b_n,state->inertia); // local angular velocity at time n
+ Quaternionr dotQ_n=DotQ(angVel_b_n,state->ori); // dQ/dt at time n
+ Quaternionr Q_half = state->ori + dt/2 * dotQ_n; // Q at time n+1/2
+ state->prevAngMom+=dt*M; // global angular momentum at time n+1/2
+ Vector3r l_b_half = state->ori.Conjugate().Rotate(state->prevAngMom); // local angular momentum at time n+1/2
+ Vector3r angVel_b_half = diagDiv(l_b_half,state->inertia); // local angular velocity at time n+1/2
+ Quaternionr dotQ_half=DotQ(angVel_b_half,Q_half); // dQ/dt at time n+1/2
+ state->ori+=dt*dotQ_half; // Q at time n+1
+ state->ori.Normalize();
+ state->angVel=state->ori.Rotate(angVel_b_half); // global angular velocity at time n+1/2
+ //if(rb->isClump()) static_cast<Clump*>(rb.get())->moveMembers();
+ static_cast<Clump*>(rb.get())->moveMembers();
+
+ LOG_TRACE("\nforce: " << F << " torque: " << M
+ << "\nglobal angular momentum at time n: " << l_n
+ << "\nlocal angular momentum at time n: " << l_b_n
+ << "\nlocal angular velocity at time n: " << angVel_b_n
+ << "\ndQ/dt at time n: " << dotQ_n
+ << "\nQ at time n+1/2: " << Q_half
+ << "\nglobal angular momentum at time n+1/2: " << state->prevAngMom
+ << "\nlocal angular momentum at time n+1/2: " << l_b_half
+ << "\nlocal angular velocity at time n+1/2: " << angVel_b_half
+ << "\ndQ/dt at time n+1/2: " << dotQ_half
+ << "\nQ at time n+1: " << state->ori
+ << "\nglobal angular velocity at time n+1/2: " << state->angVel
+ );
+}
+
+Quaternionr NewtonsDampedLaw::DotQ(const Vector3r& angVel, const Quaternionr& Q){
+ Quaternionr dotQ(Quaternionr::ZERO);
+ dotQ[0] = (-Q[1]*angVel[0]-Q[2]*angVel[1]-Q[3]*angVel[2])/2;
+ dotQ[1] = ( Q[0]*angVel[0]-Q[3]*angVel[1]+Q[2]*angVel[2])/2;
+ dotQ[2] = ( Q[3]*angVel[0]+Q[0]*angVel[1]-Q[1]*angVel[2])/2;
+ dotQ[3] = (-Q[2]*angVel[0]+Q[1]*angVel[1]+Q[0]*angVel[2])/2;
+ return dotQ;
+}
=== modified file 'pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp'
--- pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp 2009-11-23 15:36:10 +0000
+++ pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp 2009-11-26 10:34:14 +0000
@@ -41,25 +41,30 @@
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, State* clumpRBP);
bool haveBins;
+ void accurateRigidBodyRotationIntegrator(MetaBody* ncb, const shared_ptr<Body>& rb);
+ Quaternionr DotQ(const Vector3r& angVel, const Quaternionr& Q);
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;
+ /// Enable of the accurate rigid body rotation integrator
+ bool accRigidBodyRot;
#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), accRigidBodyRot(false){
#ifdef YADE_OPENMP
threadMaxVelocitySq.resize(omp_get_max_threads());
#endif
}
- REGISTER_ATTRIBUTES(StandAloneEngine,(damping)(maxVelocitySq));
+ REGISTER_ATTRIBUTES(StandAloneEngine,(damping)(maxVelocitySq)(accRigidBodyRot));
REGISTER_CLASS_AND_BASE(NewtonsDampedLaw,StandAloneEngine);
+ DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(NewtonsDampedLaw);
=== modified file 'pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp'
--- pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-11-21 16:46:58 +0000
+++ pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-11-26 10:34:14 +0000
@@ -49,8 +49,8 @@
else if(rec=="facets") recActive[REC_FACETS]=true;
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). Ignored.");
+ else if(rec=="intr") recActive[REC_INTR]=true;
+ else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: spheres, velocity, facets, colors, cpm, intr). Ignored.");
}
// cpm needs interactions
if(recActive[REC_CPM]) recActive[REC_INTR]=true;
@@ -253,6 +253,3 @@
//writer->Write();
}
-
-YADE_REQUIRE_FEATURE(PHYSPAR);
-
=== modified file 'pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp'
--- pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp 2009-10-21 15:22:14 +0000
+++ pkg/dem/Engine/StandAloneEngine/VTKRecorder.hpp 2009-11-26 10:34:14 +0000
@@ -20,7 +20,7 @@
virtual void action(MetaBody*);
private:
- REGISTER_ATTRIBUTES(PeriodicEngine,(recorders)(fileName)(compress)(skipNondynamic));
+ REGISTER_ATTRIBUTES(PeriodicEngine,(recorders)(fileName)(compress)(skipNondynamic)(skipFacetIntr));
REGISTER_CLASS_AND_BASE(VTKRecorder,PeriodicEngine);
DECLARE_LOGGER;
};
