yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2253
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Tue 2010-05-25 13:10:22 +0200
message:
  Missing files in r2249.
  
  Corresponding to : #3 - Implement homothetic=2 properly. (homothetic=1 is removed)
modified:
  pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp
  pkg/dem/Engine/GlobalEngine/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 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp'
--- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp	2010-05-04 13:56:05 +0000
+++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp	2010-05-25 11:10:22 +0000
@@ -11,9 +11,7 @@
 #include<yade/pkg-dem/Clump.hpp>
 #include<yade/pkg-common/VelocityBins.hpp>
 #include<yade/lib-base/Math.hpp>
-
-
-
+		
 YADE_PLUGIN((NewtonIntegrator));
 CREATE_LOGGER(NewtonIntegrator);
 void NewtonIntegrator::cundallDamp(const Real& dt, const Vector3r& N, const Vector3r& V, Vector3r& A){
@@ -65,8 +63,6 @@
 {
 	scene->forces.sync();
 	Real dt=scene->dt;
-	// precompute transformation increment; using Cell::getTrsfInc would get increment from the previous step, which is not right... (?) -> should be ok (B.)
-	if(scene->isPeriodic && homotheticCellResize) cellTrsfInc=dt*scene->cell->velGrad;
 	// account for motion of the periodic boundary, if we remember its last position
 	// its velocity will count as max velocity of bodies
 	// otherwise the collider might not run if only the cell were changing without any particle motion
@@ -173,27 +169,22 @@
 		FOREACH(const Real& thrMaxVSq, threadMaxVelocitySq) { maxVelocitySq=max(maxVelocitySq,thrMaxVSq); }
 	#endif
 	if(haveBins) velocityBins->binVelSqFinalize();
-	if(scene->isPeriodic) { prevCellSize=scene->cell->getSize(); }
+	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 )
+inline void NewtonIntegrator::leapfrogTranslate(Scene* scene, State* state, const body_id_t& id, const Real& dt)
 {
 	blockTranslateDOFs(state->blockedDOFs, state->accel);
 	state->vel+=dt*state->accel;
+
+	if (scene->forces.getMoveRotUsed()) state->pos+=scene->forces.getMove(id);
+	if (homotheticCellResize) {
+		// 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).
+		Vector3r dVel((scene->cell->velGrad-prevVelGrad)*/*scene->cell->wrapShearedPt(*/state->pos/*)*/); state->vel+=dVel;
+	}
 	state->pos += state->vel*dt;
-	if(scene->forces.getMoveRotUsed()) state->pos+=scene->forces.getMove(id);
-	assert(homotheticCellResize>=0 && homotheticCellResize<=2);
-	if(homotheticCellResize>0){
-		//Vector3r dPos(scene->cell->getTrsfInc()*scene->cell->wrapShearedPt(state->pos));
-		Vector3r dPos(cellTrsfInc*scene->cell->wrapShearedPt(state->pos));
-		// apply cell deformation
-		//if(homotheticCellResize>=1) it is the same test again, useless
-		state->pos+=dPos;
-		// update velocity for usage in rate dependant equations (e.g. viscous law)
-		// FIXME : it is not recommended to do that because it impacts the dynamics (this modified velocity will be used as reference in the next time-step)
-		if(homotheticCellResize==2) state->vel+=dPos/dt;
-	}
-
 }
 
 inline void NewtonIntegrator::leapfrogSphericalRotate(Scene* scene, State* state, const body_id_t& id, const Real& dt )

=== modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp'
--- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp	2010-04-25 13:18:11 +0000
+++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.hpp	2010-05-25 11:10:22 +0000
@@ -58,8 +58,6 @@
 	bool cellChanged;
 
 	public:
-		//! Store transformation increment for the current step (updated automatically)
-		Matrix3r cellTrsfInc;
 		#ifdef YADE_OPENMP
 			vector<Real> threadMaxVelocitySq;
 		#endif
@@ -70,7 +68,8 @@
 		((Real,damping,0.2,"damping coefficient for Cundall's non viscous damping (see [Chareyre2005]_) [-]"))
 		((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 aspherical bodies *only*, using formulation from [Allen1989]_, pg. 89."))
-		((int,homotheticCellResize,((void)"disabled",0),"Enable artificially moving all bodies with the periodic cell, such that its resizes are isotropic. 0: disabled, 1: position update, 2: velocity update."))
+		((bool,homotheticCellResize,false,"Enable artificially moving all bodies with the periodic cell, such that its resizes are homogeneous. The move is reflecting changes in :yref:`Cell::velGrad`, using :yref:`NewtonIntegrator::prevVelGrad`."))
+		((Matrix3r,prevVelGrad,Matrix3r::Zero(),"Store previous velocity gradient (:yref:`Cell::velGrad`) to track acceleration. |yupdate|"))
 		((vector<shared_ptr<BodyCallback> >,callbacks,,"List (std::vector in c++) of :yref:`BodyCallbacks<BodyCallback>` which will be called for each body as it is being processed."))
 		,
 		/*ctor*/