yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3696
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
timestamp: Wed 2015-07-01 20:36:34 +0200
message:
  apply the spin of the velocity gradient on particles in periodic BCs
modified:
  core/Cell.hpp
  pkg/dem/NewtonIntegrator.cpp
  pkg/dem/NewtonIntegrator.hpp


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'core/Cell.hpp'
--- core/Cell.hpp	2015-05-22 05:46:49 +0000
+++ core/Cell.hpp	2015-07-01 18:36:34 +0000
@@ -141,6 +141,7 @@
 	Matrix3r getRotation() { Matrix3r R,U; computePolarDecOfDefGrad(R,U); return R; }
 	Matrix3r getLeftStretch() { Matrix3r R,U; computePolarDecOfDefGrad(R,U); return U; }
 	Matrix3r getRightStretch() { Matrix3r R,U; computePolarDecOfDefGrad(R,U); return trsf*R.transpose(); }
+	Vector3r getSpin() {Matrix3r R=.5*(velGrad-velGrad.transpose()); return Vector3r(-R(1,2),R(0,2),-R(0,1));}
 
 	// stress measures
 	//Matrix3r getStress() { return Shop::getStress(); }
@@ -157,7 +158,7 @@
 		((Matrix3r,velGrad,Matrix3r::Zero(),,"[overridden below]"))
 		((Matrix3r,nextVelGrad,Matrix3r::Zero(),Attr::readonly,"see :yref:`Cell.velGrad`."))
 		((Matrix3r,prevVelGrad,Matrix3r::Zero(),Attr::readonly,"Velocity gradient in the previous step."))
-		((bool,homoDeform,true,,"Deform (:yref:`velGrad<Cell.velGrad>`) the cell homothetically, by adjusting positions and velocities of bodies. The velocity change is obtained by deriving the expression v=∇v.x, where ∇v is the macroscopic velocity gradient, giving in an incremental form: Δv=Δ ∇v x + ∇v Δx. As a result, velocities are modified as soon as ``velGrad`` changes, according to the first term: Δv(t)=Δ ∇v x(t), while the 2nd term reflects a convective term: Δv'= ∇v v(t-dt/2)."))
+		((int,homoDeform,2,,"If >0, deform (:yref:`velGrad<Cell.velGrad>`) the cell homothetically by adjusting positions and velocities of bodies. The velocity change is obtained by deriving the expression v=∇v.x, where ∇v is the macroscopic velocity gradient, giving in an incremental form: Δv=Δ ∇v x + ∇v Δx. As a result, velocities are modified as soon as ``velGrad`` changes, according to the first term: Δv(t)=Δ ∇v x(t), while the 2nd term reflects a convective term: Δv'= ∇v v(t-dt/2). The second term is neglected if homoDeform=1. All terms are included if homoDeform=2 (default)"))
 		((bool,velGradChanged,false,Attr::readonly,"true when velGrad has been changed manually (see also :yref:`Cell.nextVelGrad`)")),
 		/*ctor*/ _invTrsf=Matrix3r::Identity(); integrateAndUpdate(0),
 		/*py*/
@@ -188,6 +189,7 @@
 		.def("getRotation",&Cell::getRotation,"Returns rotation of the cell (orthogonal matrix $\\mat{R}$ from polar decomposition $\\mat{F}=\\mat{RU}$ )")
 		.def("getLeftStretch",&Cell::getLeftStretch,"Returns left (spatial) stretch tensor of the cell (matrix $\\mat{U}$ from polar decomposition $\\mat{F}=\\mat{RU}$ )")
 		.def("getRightStretch",&Cell::getRightStretch,"Returns right (material) stretch tensor of the cell (matrix $\\mat{V}$ from polar decomposition $\\mat{F}=\\mat{RU}=\\mat{VR}\\ \\rightarrow\\ \\mat{V}=\\mat{FR}^T$ )")
+		.def("getSpin",&Cell::getSpin,"Returns the spin defined by the skew symmetric part of :yref:`velGrad<Cell.velGrad>`") 
 		.def_readonly("shearTrsf",&Cell::_shearTrsf,"Current skew+rot transformation (no resize)")
 		.def_readonly("unshearTrsf",&Cell::_unshearTrsf,"Inverse of the current skew+rot transformation (no resize)")
 		.add_property("hSize0",&Cell::getHSize0,"Value of untransformed hSize, with respect to current :yref:`trsf<Cell::trsf>` (computed as :yref:`trsf<Cell::trsf>` ⁻¹ × :yref:`hSize<Cell::hSize>`.")

=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp	2015-04-24 15:54:07 +0000
+++ pkg/dem/NewtonIntegrator.cpp	2015-07-01 18:36:34 +0000
@@ -96,6 +96,8 @@
 		scene->cell->velGradChanged=0; scene->cell->nextVelGrad=Matrix3r::Zero();}
 	homoDeform=scene->cell->homoDeform;
 	dVelGrad=scene->cell->velGrad-prevVelGrad;
+	Matrix3r R=.5*(dVelGrad-dVelGrad.transpose());
+	dSpin=Vector3r(-R(1,2),R(0,2),-R(0,1));
 	// 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
@@ -165,7 +167,7 @@
 				if (densityScaling) linAccel*=state->densityScaling;
 				if(state->isDamped) cundallDamp2nd(dt,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 (isPeriodic && homoDeform) linAccel+=prevVelGrad*state->vel;
+				if (isPeriodic && homoDeform>1) linAccel+=prevVelGrad*state->vel;
 				//finally update velocity
 				state->vel+=dt*linAccel;
 				// angular acceleration
@@ -179,7 +181,7 @@
 					if(state->isDamped) cundallDamp1st(m,state->angVel);
 				}
 			// reflect macro-deformation even for non-dynamic bodies
-			} else if (isPeriodic && homoDeform) state->vel+=dt*prevVelGrad*state->vel;
+			} else if (isPeriodic && homoDeform>1) state->vel+=dt*prevVelGrad*state->vel;
 
 			// update positions from velocities (or torque, for the aspherical integrator)
 			leapfrogTranslate(state,id,dt);
@@ -219,6 +221,7 @@
 
 void NewtonIntegrator::leapfrogSphericalRotate(State* state, const Body::id_t& id, const Real& dt )
 {
+ 	if(scene->isPeriodic && homoDeform) {state->angVel+=dSpin;}
 	Real angle2=state->angVel.squaredNorm();
 	if (angle2!=0 and ( (mask<=0) or ((mask>0) and (Body::byId(id)->maskCompatible(mask))) )) {//If we have an angular velocity, we make a rotation
 		Real angle=sqrt(angle2);
@@ -236,6 +239,8 @@
 }
 
 void NewtonIntegrator::leapfrogAsphericalRotate(State* state, const Body::id_t& id, const Real& dt, const Vector3r& M){
+	//FIXME: where to increment angular velocity like this? Only done for spherical rotations at the moment
+	//if(scene->isPeriodic && homoDeform) {state->angVel+=dSpin;}
 	Matrix3r A=state->ori.conjugate().toRotationMatrix(); // rotation matrix from global to local r.f.
 	const Vector3r l_n = state->angMom + dt/2. * M; // global angular momentum at time n
 	const Vector3r l_b_n = A*l_n; // local angular momentum at time n

=== modified file 'pkg/dem/NewtonIntegrator.hpp'
--- pkg/dem/NewtonIntegrator.hpp	2015-03-05 08:32:27 +0000
+++ pkg/dem/NewtonIntegrator.hpp	2015-07-01 18:36:34 +0000
@@ -40,11 +40,12 @@
 	#endif
 	// whether the cell has changed from the previous step
 	bool cellChanged;
-	bool homoDeform;
+	int homoDeform;
 	
 	// wether a body has been selected in Qt view
 	bool bodySelected;
 	Matrix3r dVelGrad;
+	Vector3r dSpin;
 
 	public:
 		bool densityScaling;// internal for density scaling