yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2922
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Tue 2011-09-20 12:58:18 +0200
message:
  - fix the "homoDeform" bug https://bugs.launchpad.net/yade/+bug/836867, and adapt regression tests so they don't test homoDeform=1
  - add gravity to Newton
  - tune front page of the documentation to be more google-scholar-friendly
modified:
  core/Cell.cpp
  core/Cell.hpp
  doc/sphinx/conf.py
  pkg/dem/NewtonIntegrator.cpp
  pkg/dem/NewtonIntegrator.hpp
  py/tests/pbc.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 'core/Cell.cpp'
--- core/Cell.cpp	2011-01-30 15:43:18 +0000
+++ core/Cell.cpp	2011-09-20 10:58:18 +0000
@@ -32,8 +32,6 @@
 	_hasShear=(hSize(0,1)!=0 || hSize(0,2)!=0 || hSize(1,0)!=0 || hSize(1,2)!=0 || hSize(2,0)!=0 || hSize(2,1)!=0);
 	// OpenGL shear matrix (used frequently)
 	fillGlShearTrsfMatrix(_glShearTrsfMatrix);
-
-	if(!(homoDeform==HOMO_NONE || homoDeform==HOMO_POS || homoDeform==HOMO_VEL || homoDeform==HOMO_VEL_2ND)) throw std::invalid_argument("Cell.homoDeform must be in {0,1,2,3}.");
 }
 
 void Cell::fillGlShearTrsfMatrix(double m[16]){

=== modified file 'core/Cell.hpp'
--- core/Cell.hpp	2011-09-02 19:23:13 +0000
+++ core/Cell.hpp	2011-09-20 10:58:18 +0000
@@ -87,9 +87,9 @@
 
 	// relative position and velocity for interaction accross multiple cells
 	Vector3r intrShiftPos(const Vector3i& cellDist) const { return hSize*cellDist.cast<Real>(); }
-	Vector3r intrShiftVel(const Vector3i& cellDist) const { if(homoDeform==HOMO_VEL || homoDeform==HOMO_VEL_2ND) return velGrad*hSize*cellDist.cast<Real>(); return Vector3r::Zero(); }
+	Vector3r intrShiftVel(const Vector3i& cellDist) const { return velGrad*hSize*cellDist.cast<Real>(); }
 	// return body velocity while taking away mean field velocity (coming from velGrad) if the mean field velocity is applied on velocity
-	Vector3r bodyFluctuationVel(const Vector3r& pos, const Vector3r& vel, const Matrix3r& prevVelGrad) const { if(homoDeform==HOMO_VEL || homoDeform==HOMO_VEL_2ND) return (vel-prevVelGrad*pos); return vel; }
+	Vector3r bodyFluctuationVel(const Vector3r& pos, const Vector3r& vel, const Matrix3r& prevVelGrad) const { return (vel-prevVelGrad*pos); }
 
 	// get/set current shape; setting resets trsf to identity
 	Matrix3r getHSize() const { return hSize; }
@@ -133,8 +133,10 @@
 		/* normal attributes */
 		((Matrix3r,velGrad,Matrix3r::Zero(),,"Velocity gradient of the transformation; used in :yref:`NewtonIntegrator`. Values of :yref:`velGrad<Cell.velGrad>` accumulate in :yref:`trsf<Cell.trsf>` at every step."))
 		((Matrix3r,prevVelGrad,Matrix3r::Zero(),Attr::readonly,"Velocity gradient in the previous step."))
-		((int,homoDeform,3,Attr::triggerPostLoad,"Deform (:yref:`velGrad<Cell.velGrad>`) the cell homothetically, by adjusting positions or velocities of particles. The values have the following meaning: 0: no homothetic deformation, 1: set absolute particle positions directly (when ``velGrad`` is non-zero), but without changing their velocity, 2: adjust particle velocity (only when ``velGrad`` changed) with Δv_i=Δ ∇v x_i. 3: as 2, but include a 2nd order term in addition -- the derivative of 1 (convective term in the velocity update).")),
-		/*deprec*/ ((Hsize,hSize,"conform to Yade's names convention.")),
+		((int,homoDeformOld,3,,"[DEPRECATED: homoDeform=3 is the only option left, kept here for compatibility] Deform (:yref:`velGrad<Cell.velGrad>`) the cell homothetically, by adjusting positions or velocities of particles. The values have the following meaning: 0: no homothetic deformation, 1: set absolute particle positions directly (when ``velGrad`` is non-zero), but without changing their velocity, 2: adjust particle velocity (only when ``velGrad`` changed) with Δv_i=Δ ∇v x_i. 3: as 2, but include a 2nd order term in addition -- the derivative of 1 (convective term in the velocity update).")),
+		/*deprec*/
+		((Hsize,hSize,"conform to Yade's names convention."))
+		((homoDeform,homoDeformOld,"[DEPRECATED: homoDeform=3 is the only option left, kept here for compatibility]")),
 		/*init*/ ,
 		/*ctor*/ _invTrsf=Matrix3r::Identity(); integrateAndUpdate(0),
 		/*py*/

=== modified file 'doc/sphinx/conf.py'
--- doc/sphinx/conf.py	2011-06-08 16:21:31 +0000
+++ doc/sphinx/conf.py	2011-09-20 10:58:18 +0000
@@ -648,9 +648,6 @@
 \vspace{20 mm}
 \text{\sffamily\bfseries\Huge Yade Documentation}\\
 \vspace{5 mm}
-\textit{\sffamily\huge Release '''\
-+yade.config.revision\
-+r'''}\\
 \vspace{70 mm}
 \begin{sffamily}\bfseries\Large
 Václav Šmilauer, Emanuele Catalano, Bruno Chareyre, Sergei Dorofeenko, Jerome Duriez, Anton Gladky, Janek Kozicki, Chiara Modenese, Luc Scholtès, Luc Sibille, Jan Stránský, Klaus Thoeni
@@ -660,7 +657,9 @@
 
 \vfill
 % Bottom of the page
-{\large \today}
+\textit{\Large Release '''\
++yade.config.revision\
++r''', \today}
 \end{flushright}
 
 \end{titlepage}

=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp	2011-09-02 19:23:13 +0000
+++ pkg/dem/NewtonIntegrator.cpp	2011-09-20 10:58:18 +0000
@@ -73,7 +73,7 @@
 	bodySelected=(scene->selectedBody>=0);
 	if(warnNoForceReset && scene->forces.lastReset<scene->iter) LOG_WARN("O.forces last reset in step "<<scene->forces.lastReset<<", while the current step is "<<scene->iter<<". Did you forget to include ForceResetter in O.engines?");
 	const Real& dt=scene->dt;
-	homoDeform=(scene->isPeriodic ? scene->cell->homoDeform : -1); // -1 for aperiodic simulations
+
 	dVelGrad=scene->cell->velGrad-prevVelGrad;
 	// account for motion of the periodic boundary, if we remember its last position
 	// its velocity will count as max velocity of bodies
@@ -121,8 +121,8 @@
 				scene->forces.addForce(id, f);
 				#endif
 			}
-			//in most cases, the initial force on clumps will zero and next line is not changing f and m, but make sure we don't miss something (e.g. user defined forces on clumps)
-			f=scene->forces.getForce(id); m=scene->forces.getTorque(id);
+			//in most cases, the initial force on clumps will be zero and next line is not changing f and m, but make sure we don't miss something (e.g. user defined forces on clumps)
+			f=scene->forces.getForce(id)+gravity; m=scene->forces.getTorque(id);
 			#ifdef YADE_DEBUG
 				if(isnan(f[0])||isnan(f[1])||isnan(f[2])) throw runtime_error(("NewtonIntegrator: NaN force acting on #"+lexical_cast<string>(id)+".").c_str());
 				if(isnan(m[0])||isnan(m[1])||isnan(m[2])) throw runtime_error(("NewtonIntegrator: NaN torque acting on #"+lexical_cast<string>(id)+".").c_str());
@@ -148,7 +148,7 @@
 				Vector3r linAccel=computeAccel(f,state->mass,state->blockedDOFs);
 				cundallDamp2nd(dt,f,fluctVel,linAccel);
 				//This is the convective term, appearing in the time derivation of Cundall/Thornton expression (dx/dt=velGrad*pos -> d²x/dt²=dvelGrad/dt*pos+velGrad*vel), negligible in many cases but not for high speed large deformations (gaz or turbulent flow).
-				if (homoDeform==Cell::HOMO_VEL_2ND) linAccel+=prevVelGrad*state->vel;
+				linAccel+=prevVelGrad*state->vel;
 				//finally update velocity
 				state->vel+=dt*linAccel;
 				// angular acceleration
@@ -188,16 +188,12 @@
 
 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.
-		//NOTE : if the velocity is updated before moving the body, it means the current velGrad (i.e. before integration in cell->integrateAndUpdate) will be effective for the current time-step. Is it correct? If not, this velocity update can be moved just after "state->pos += state->vel*dt", meaning the current velocity impulse will be applied at next iteration, after the contact law. (All this assuming the ordering is resetForces->integrateAndUpdate->contactLaw->PeriCompressor->NewtonsLaw. Any other might fool us.)
-		//NOTE : dVel defined without wraping the coordinates means bodies out of the (0,0,0) period can move realy fast. It has to be compensated properly in the definition of relative velocities (see Ig2 functors and contact laws).
+	// update velocity reflecting changes in the macroscopic velocity field, making the problem homothetic.
+	//NOTE : if the velocity is updated before moving the body, it means the current velGrad (i.e. before integration in cell->integrateAndUpdate) will be effective for the current time-step. Is it correct? If not, this velocity update can be moved just after "state->pos += state->vel*dt", meaning the current velocity impulse will be applied at next iteration, after the contact law. (All this assuming the ordering is resetForces->integrateAndUpdate->contactLaw->PeriCompressor->NewtonsLaw. Any other might fool us.)
+	//NOTE : dVel defined without wraping the coordinates means bodies out of the (0,0,0) period can move realy fast. It has to be compensated properly in the definition of relative velocities (see Ig2 functors and contact laws).
 		//Reflect mean-field (periodic cell) acceleration in the velocity
-		Vector3r dVel=dVelGrad*state->pos;
-		state->vel+=dVel;
-	} else if (homoDeform==Cell::HOMO_POS){
-		state->pos+=scene->cell->velGrad*state->pos*dt;
-	}
+	Vector3r dVel=dVelGrad*state->pos;
+	state->vel+=dVel;
 	if (!bodySelected || scene->selectedBody!=id) state->pos+=state->vel*dt;
 }
 

=== modified file 'pkg/dem/NewtonIntegrator.hpp'
--- pkg/dem/NewtonIntegrator.hpp	2011-04-22 09:01:59 +0000
+++ pkg/dem/NewtonIntegrator.hpp	2011-09-20 10:58:18 +0000
@@ -40,8 +40,7 @@
 	bool cellChanged;
 	// wether a body has been selected in Qt view
 	bool bodySelected;
-	int homoDeform; // updated from scene at every call; -1 for aperiodic simulations, otherwise equal to scene->cell->homoDeform
-	Matrix3r dVelGrad; // dtto
+	Matrix3r dVelGrad;
 
 	public:
 		// function to save maximum velocity, for the verlet-distance optimization
@@ -54,6 +53,7 @@
 		virtual void action();
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR(NewtonIntegrator,GlobalEngine,"Engine integrating newtonian motion equations.",
 		((Real,damping,0.2,,"damping coefficient for Cundall's non viscous damping (see [Chareyre2005]_) [-]"))
+		((Vector3r,gravity,Vector3r::Zero(),,"Gravitational acceleration (effectifely replaces GravityEngine)."))
 		((Real,maxVelocitySq,NaN,,"store square of max. velocity, for informative purposes; computed again at every step. |yupdate|"))
 		((bool,exactAsphericalRot,true,,"Enable more exact body rotation integrator for :yref:`aspherical bodies<Body.aspherical>` *only*, using formulation from [Allen1989]_, pg. 89."))
 		((Matrix3r,prevVelGrad,Matrix3r::Zero(),,"Store previous velocity gradient (:yref:`Cell::velGrad`) to track acceleration. |yupdate|"))

=== modified file 'py/tests/pbc.py'
--- py/tests/pbc.py	2011-01-31 19:09:24 +0000
+++ py/tests/pbc.py	2011-09-20 10:58:18 +0000
@@ -29,7 +29,6 @@
 		O.bodies[1].state.vel=self.initVel
 		O.engines=[NewtonIntegrator(warnNoForceReset=False)]
 		O.cell.velGrad=Matrix3(0,0,0, 0,0,0, 0,0,-1)
-		O.cell.homoDeform=3
 		O.dt=0 # do not change positions with dt=0 in NewtonIntegrator, but still update velocities from velGrad
 	def testVelGrad(self):
 		'PBC: velGrad changes hSize but not hSize0, accumulates in trsf'
@@ -58,29 +57,13 @@
 		O.cell.setBox(2.55,11,45)
 		self.assert_(O.cell.hSize==Matrix3(2.55,0,0, 0,11,0, 0,0,45));
 	def testHomotheticResizeVel(self):
-		"PBC: homothetic cell deformation adjusts particle velocity (homoDeform==3)"
+		"PBC: homothetic cell deformation adjusts particle velocity "
 		O.dt=1e-5
 		O.step()
 		s1=O.bodies[1].state
 		self.assertAlmostEqual(s1.vel[2],self.initVel[2]+self.initPos[2]*O.cell.velGrad[2,2])
-	def testHomotheticResizePos(self):
-		"PBC: homothetic cell deformation adjusts particle position (homoDeform==1)"
-		O.cell.homoDeform=1
-		O.step()
-		s1=O.bodies[1].state
-		self.assertAlmostEqual(s1.vel[2],self.initVel[2])
-		self.assertAlmostEqual(s1.pos[2],self.initPos[2]+self.initPos[2]*O.cell.velGrad[2,2]*O.dt)
 	def testScGeomIncidentVelocity(self):
-		"PBC: ScGeom computes incident velocity correctly (homoDeform==3)"
-		O.step()
-		O.engines=[InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[])]
-		i=utils.createInteraction(0,1)
-		self.assertEqual(self.initVel,i.geom.incidentVel(i,avoidGranularRatcheting=True))
-		self.assertEqual(self.initVel,i.geom.incidentVel(i,avoidGranularRatcheting=False))
-		self.assertAlmostEqual(self.relDist[1],1-i.geom.penetrationDepth)
-	def testScGeomIncidentVelocity_homoPos(self):
-		"PBC: ScGeom computes incident velocity correctly (homoDeform==1)"
-		O.cell.homoDeform=1
+		"PBC: ScGeom computes incident velocity correctly"
 		O.step()
 		O.engines=[InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[])]
 		i=utils.createInteraction(0,1)
@@ -88,31 +71,20 @@
 		self.assertEqual(self.initVel,i.geom.incidentVel(i,avoidGranularRatcheting=False))
 		self.assertAlmostEqual(self.relDist[1],1-i.geom.penetrationDepth)
 	def testL3GeomIncidentVelocity(self):
-		"PBC: L3Geom computes incident velocity correctly (homoDeform==3)"
+		"PBC: L3Geom computes incident velocity correctly"
 		O.step()
 		O.engines=[ForceResetter(),InteractionLoop([Ig2_Sphere_Sphere_L3Geom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L3Geom_FrictPhys_ElPerfPl(noBreak=True)]),NewtonIntegrator()]
 		i=utils.createInteraction(0,1) 
 		O.dt=1e-10; O.step() # tiny timestep, to not move the normal too much
 		self.assertAlmostEqual(self.initVel.norm(),(i.geom.u/O.dt).norm())
-	def testL3GeomIncidentVelocity_homoPos(self):
-		"PBC: L3Geom computes incident velocity correctly (homoDeform==1)"
-		O.cell.homoDeform=1; O.step()
-		O.engines=[ForceResetter(),InteractionLoop([Ig2_Sphere_Sphere_L3Geom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_L3Geom_FrictPhys_ElPerfPl(noBreak=True)]),NewtonIntegrator()]
-		i=utils.createInteraction(0,1) 
-		O.dt=1e-10; O.step()
-		self.assertAlmostEqual(self.initVel.norm(),(i.geom.u/O.dt).norm())
-		#self.assertAlmostEqual(self.relDist[1],1-i.geom.penetrationDepth)
 	def testKineticEnergy(self):
-		"PBC: utils.kineticEnergy considers only fluctuation velocity, not the velocity gradient (homoDeform==3)"
+		"PBC: utils.kineticEnergy considers only fluctuation velocity, not the velocity gradient"
 		O.step() # updates velocity with homotheticCellResize
 		# ½(mv²+ωIω)
 		# #0 is still, no need to add it; #1 has zero angular velocity
 		# we must take self.initVel since O.bodies[1].state.vel now contains the homothetic resize which utils.kineticEnergy is supposed to compensate back 
 		Ek=.5*O.bodies[1].state.mass*self.initVel.squaredNorm()
 		self.assertAlmostEqual(Ek,utils.kineticEnergy())
-	def testKineticEnergy_homoPos(self):
-		"PBC: utils.kineticEnergy considers only fluctuation velocity, not the velocity gradient (homoDeform==1)"
-		O.cell.homoDeform=1; O.step()
-		self.assertAlmostEqual(.5*O.bodies[1].state.mass*self.initVel.squaredNorm(),utils.kineticEnergy())
+