yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2637
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Mon 2011-01-03 14:24:40 +0100
message:
  - Don't damp clump members. Instead, compute undamped acceleration and damp it at the clump's level.
  - Similar mechanism implemented for DOF's blocking : compute unblocked accelerations and block at the end (remember blocking one DOF triggers tests and divisions).
  - Still possible optimizations (especially for spherical rotations and translations), see https://blueprints.launchpad.net/yade/+spec/clump-forces.
  - One FIXME : aspherical damping is not correct a.t.m. since it uses torque where angAccel should be provided (suggested fix in sources).
  - Probably fixes a few inconsistencies in previous commits.
modified:
  pkg/dem/NewtonIntegrator.cpp


--
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/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp	2011-01-03 11:22:51 +0000
+++ pkg/dem/NewtonIntegrator.cpp	2011-01-03 13:24:40 +0000
@@ -19,6 +19,12 @@
 	for(int i=0; i<3; i++) A[i]*= 1 - damping*Mathr::Sign ( N[i]*(V[i] + 0.5*dt*A[i]) );
 }
 
+void NewtonIntegrator::blockDOFs(State* state, const bool& rotational){
+	if(unlikely(state->blockedDOFs!=0))
+	for(int i=0; i<3; i++) if(state->blockedDOFs & State::axisDOF(i,false)){
+		if (rotational) state->angAccel[i]=0; else state->accel[i]=0;}
+}
+
 Vector3r NewtonIntegrator::computeAccel(const Vector3r& force, const Real& mass, int blockedDOFs){
 	if(likely(blockedDOFs==0)) return force/mass;
 	Vector3r ret(Vector3r::Zero());
@@ -34,22 +40,23 @@
 
 Vector3r NewtonIntegrator::clumpMemberAccel(const Body::id_t& memberId, const State* memberState, const State* clumpState){
 	const Vector3r& f=scene->forces.getForce(memberId);
-	Vector3r dA=computeAccel(f,clumpState->mass,clumpState->blockedDOFs);
-	cundallDamp(scene->dt,f,memberState->vel,dA); // damp using velocities of the member
+	Vector3r dA=computeAccel(f,clumpState->mass,0);//don't block yet, use the fast version and block the final clump acceleration
+// 	cundallDamp(scene->dt,f,memberState->vel,dA); // damp using velocities of the member
 	return dA;
 }
 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 dAA=computeAngAccel(m+(memberState->pos-clumpState->pos).cross(f),clumpState->inertia,clumpState->blockedDOFs);
-	cundallDamp(scene->dt,m,memberState->angVel,dAA);
+	Vector3r dAA=computeAngAccel(m+(memberState->pos-clumpState->pos).cross(f),clumpState->inertia,0);
+// 	cundallDamp(scene->dt,m,memberState->angVel,dAA);
 	return dAA;
 }
 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);
 	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
+// 	for(int i=0; i<3; i++) if(clumpState->blockedDOFs & State::axisDOF(i,true)) dM[i]=0.;
+	//FIXME(a) see FIXME(b)
+// 	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;
 }
 
@@ -150,6 +157,7 @@
 					leapfrogSphericalRotate(state,id,dt);
 				} else { // exactAsphericalRot enabled & aspherical body
 					// damp the torque
+					//FIXME(b) : cundallDamp will be multiplying mm by dt, as if it were an acceleration! Divide by mass here? Another similar occurence is below for clumps (FIXME(c))
 					Vector3r mm(m); cundallDamp(0,m,state->angVel,mm);
 					leapfrogAsphericalRotate(state,id,dt,mm);
 				}
@@ -157,7 +165,9 @@
 				// 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
-				if(state->blockedDOFs!=State::DOF_ALL){ // if the clump has all DoFs blocked, forces on members would be computed uselessly
+				if(state->blockedDOFs!=State::DOF_ALL){ // if the clump has all DoFs blocked, forces on members and accelerations would be computed uselessly
+					state->accel=computeAccel(f,state->mass,0);
+					state->angAccel=computeAngAccel(M,state->inertia,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());
@@ -166,17 +176,25 @@
 						if(aspherical) M+=clumpMemberTorque(memberId,memberState,state);
 						else state->angAccel+=clumpMemberAngAccel(memberId,memberState,state);
 					}
-					state->accel=computeAccel(f,state->mass,state->blockedDOFs); cundallDamp(dt,f,fluctVel,state->accel);
+					//Translational accel.
+					blockDOFs (state,false);
+					cundallDamp(dt,state->accel,fluctVel,state->accel);
+					state->vel+=dt*state->accel;
+					//Rotational accel.
 					if(!aspherical){
-						state->angAccel+=computeAngAccel(M,state->inertia,state->blockedDOFs);
-						cundallDamp(dt,m,state->angVel,state->angAccel);
+						blockDOFs(state,true);
+						cundallDamp(dt,state->angAccel,state->angVel,state->angAccel);
 						state->angVel+=dt*state->angAccel;}
-					state->vel+=dt*state->accel;
+					//FIXME(c) this equation is replacing FIXME(a), dragging the same problem
+					else cundallDamp(dt,M,state->angVel,M);
 				}
 				// translation
 				leapfrogTranslate(state,id,dt);
 				// rotation (asphericalRotate will not keep velocity constant if State::DOF_ALL, so we skip it in this special case)
-				if(aspherical && state->blockedDOFs!=State::DOF_ALL) leapfrogAsphericalRotate(state,id,dt,M);
+				if(aspherical && state->blockedDOFs!=State::DOF_ALL) {
+					leapfrogAsphericalRotate(state,id,dt,M);
+					//blockDOFs (state,true);//proposal for handling arbitrary number of blocked DOFs?
+				}
 				else leapfrogSphericalRotate(state,id,dt);
 				// move individual members of the clump, save maxima velocity (for collider stride)
 				FOREACH(Clump::MemberMap::value_type mm, b->shape->cast<Clump>().members){