yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2629
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Mon 2010-12-27 11:20:51 +0100
message:
  -Fix kinetic energy definition.
modified:
  pkg/dem/NewtonIntegrator.cpp
  pkg/dem/Shop.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	2010-12-23 14:01:33 +0000
+++ pkg/dem/NewtonIntegrator.cpp	2010-12-27 10:20:51 +0000
@@ -130,7 +130,7 @@
 					Matrix3r mI; mI<<state->inertia[0],0,0, 0,state->inertia[1],0, 0,0,state->inertia[2];
 					Matrix3r T(state->ori);
 					Erot=.5*b->state->angVel.transpose().dot((T.transpose()*mI*T)*b->state->angVel);
-				} else { Erot=state->angVel.dot(state->inertia.cwise()*state->angVel); }
+				} else { Erot=0.5*state->angVel.dot(state->inertia.cwise()*state->angVel); }
 				if(!kinSplit) scene->energy->add(Etrans+Erot,"kinetic",kinEnergyIx,/*non-incremental*/true);
 				else{ scene->energy->add(Etrans,"kinTrans",kinEnergyTransIx,true); scene->energy->add(Erot,"kinRot",kinEnergyRotIx,true); }
 			}
@@ -146,7 +146,7 @@
 					leapfrogSphericalRotate(state,id,dt);
 				} else { // exactAsphericalRot enabled & aspherical body
 					// damp the torque
-					Vector3r mm(m); cundallDamp(0,m,state->angVel,mm); 
+					Vector3r mm(m); cundallDamp(0,m,state->angVel,mm);
 					leapfrogAsphericalRotate(state,id,dt,mm);
 				}
 			} else if (b->isClump()){

=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp	2010-12-26 15:42:43 +0000
+++ pkg/dem/Shop.cpp	2010-12-27 10:20:51 +0000
@@ -99,7 +99,7 @@
 
 	// remove all potential interactions
 	scene->interactions->eraseNonReal();
-	// adjust Interaction::cellDist for real interactions; 
+	// adjust Interaction::cellDist for real interactions;
 	FOREACH(const shared_ptr<Interaction>& i, *scene->interactions){
 		Body::id_t id1=i->getId1(),id2=i->getId2();
 		// this must be the same for both old and new interaction: cell2-cell1+cellDist
@@ -194,7 +194,7 @@
 			Matrix3r mI; mI<<state->inertia[0],0,0, 0,state->inertia[1],0, 0,0,state->inertia[2];
 			E+=.5*state->angVel.transpose().dot((T.transpose()*mI*T)*state->angVel);
 		}
-		else { E+=state->angVel.dot(state->inertia.cwise()*state->angVel);}
+		else { E+=0.5*state->angVel.dot(state->inertia.cwise()*state->angVel);}
 		if(maxId && E>maxE) { *maxId=b->getId(); maxE=E; }
 		ret+=E;
 	}
@@ -341,7 +341,7 @@
 }
 /* Project 3d point into 2d using spiral projection along given axis;
  * the returned tuple is
- * 	
+ *
  *  (height relative to the spiral, distance from axis, theta )
  *
  * dH_dTheta is the inclination of the spiral (height increase per radian),
@@ -450,10 +450,10 @@
 		if(!phys) continue;
 		if(!geom) continue;
 		const Body::id_t id1=I->getId1(), id2=I->getId2();
-		
+
 		Real minRad=(geom->refR1<=0?geom->refR2:(geom->refR2<=0?geom->refR1:min(geom->refR1,geom->refR2)));
 		Real crossSection=Mathr::PI*pow(minRad,2);
-		
+
 		Vector3r normalStress=((1./crossSection)*geom->normal.dot(phys->normalForce))*geom->normal;
 		Vector3r shearStress;
 		for(int i=0; i<3; i++){
@@ -461,10 +461,10 @@
 			shearStress[i]=geom->normal[ix1]*phys->shearForce[ix1]+geom->normal[ix2]*phys->shearForce[ix2];
 			shearStress[i]/=crossSection;
 		}
-		
+
 		bodyStates[id1].normStress+=normalStress;
 		bodyStates[id2].normStress+=normalStress;
-		
+
 		bodyStates[id1].shearStress+=shearStress;
 		bodyStates[id2].shearStress+=shearStress;
 	}
@@ -528,14 +528,14 @@
 
 		/*Real fT0=fT[0]; Real fT1=fT[1]; Real fT2=fT[2];
 		Real n0=n[0]; Real n1=n[1]; Real n2=n[2];
-		
+
 		Real s00 = l*(n0*n0*fN + fT0*n0);
 		Real s01 = l*(n0*n1*fN + .5*(fT0*n1 + fT1*n0));
 		Real s02 = l*(n0*n2*fN + .5*(fT0*n2 + fT2*n0));
 		Real s11 = l*(n1*n1*fN + fT1*n1);
 		Real s12 = l*(n1*n2*fN + .5*(fT1*n2 + fT2*n1));
 		Real s22 = l*(n2*n2*fN + fT2*n2);
-		
+
 		stress(0,0) += s00;
 		stress(0,1) += s01;
 		stress(1,0) += s01;